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Abstract 

This  thesis  develops  two  methods  of  compressing  a  track  of  radar  observations  of  a  satellite  into  a 
single  state  vector  and  associated  covariance  matrix,  and  a  method  of  estimating  orbits  using  results  from 
multiple  tracks.  The  track  compression  uses  least  squares  diSerential  correction  to  determine  a  state 
vector  at  the  central  observation  time.  One  method  uses  a  Taylor  series  expansion  in  time,  modeling  first 
and  second  order  two-bo<fy  effects  and  first  order  J2  zonal  harmonic  effects  to  represent  the  system 
dynamics.  The  second  method  uses  numerical  integration  of  the  equations  of  motion  using  two-bodjy  and 
J2  effects  to  represent  the  iq^stem  (Ramies. 

The  resulting  state  vectors  and  covariance  matrices  are  then  used  to  estimate  the  satellite’s  oibit, 
also  using  least  squares  differential  correction.  Numerical  integration  using  two-lxxfy,  J2  and  an 
atmospheric  drag  model  is  used  to  represent  the  (fynamics.  This  orbit  estimation  produces  a  state  vector 
which  includes  the  ballistic  coefficient,  as  well  as  an  associated  covariance  matrix. 

Finally,  a  one-fiftieth  scale  demonstration  of  the  fiill  AFSPC  catalog  of  satellites  and  debris  is 
conducted  to  demonstrate  the  improvement  in  accuracy  over  current  practice  which  results.  The  truth 
model  includes  J2  zonal  harmonic  effects  and  an  atmospheric  drag  model.  This  demonstration  shows  that 
the  orbits  of  90%  of  the  entire  catalog  of  objects  can  be  estimated  with  sufficient  accuracy  to  allow 
position  determination  within  one  kilometer  after  only  two  days  of  tracking.  Within  four  days,  most 
satellite  positions  are  determined  within  fiffy  meters. 

This  improvement  in  catalog  accuracy  would  result  in  a  reduced  routine  radar  tracking  workload, 
fewer  collision-avoidance  maneuvers  for  manned  space  flints,  more  accurate  over-flight  predictions, 
faster  acquisition  for  ground-based  observation  or  targeting  of  satellites  and  smaller  rendezvous  fuel 
budgets,  to  name  just  a  few  advantages. 
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ORBIT  ESTIMATION  USING  TRACK  COMPRESSION  AND 


LEAST  SQUARES  DIFFERENTIAL  CORRECTION 


/.  Introduction 


1.1  Background  [4] 

Air  Force  Space  Command  (AFSPC)  is  tasked  with  the  mission  of  tracking  all  satellites  and 
debris  in  oibit  about  the  Earth.  This  mission  is  executed  1^  the  Space  Control  Center  (SCC),  located  in 
Ch^enne  Mountain,  which  operates  the  Space  Surveillance  Networic  (SSN).  The  SSN  is  a  world-wide 
network  of  radar  sites  and  optical  tracking  cameras  which  track  satellites  as  th^  pass  overhead. 
(Throughout  this  thesis,  the  term  satellite  is  used  in  its  broader  sense,  to  denote  any  otgect  orbiting  the 
Earth,  whether  debris  or  functioning  equipment.)  The  sites  then  send  their  information  to  the  SCC  where 
it  is  centrally-processed  using  least  squares  algorithms  to  estimate  the  satellites’  orbits. 

The  radar  sites  of  the  SSN  have  the  capability  of  producing  50  to  100  observations  of  a  satellite 
per  second.  Thus  a  single  track  of  data,  firom  rise  to  set,  can  consist  of  over  100,000  observations  -  for  a 
single  pass  of  a  single  satellite  over  a  single  site.  When  this  is  factored  by  the  number  of  objects  in  the 
catalog,  and  the  number  of  sites  in  the  SSN,  the  enormous  amount  of  data  required  for  this  mission 
becomes  evident. 

1.2  Problem  Statement 

Currently,  the  SCC  does  not  possess  the  computing  resources  to  process  all  of  this  data.  Current 
practice  consists  of  the  radar  sites  reducing  the  size  of  their  observation  tracks  by  throwing  away  all  but 
three  observations.  These  three  observations,  one  near  the  beginning  of  the  track,  one  at  the  middle,  and 
one  near  the  end,  are  then  sent  to  the  SCC  for  processing.  This  reduced  workload  allows  the  SCC  to  keep 
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up  with  the  incoming  data  at  the  expense  of  estimating  orbits  with  far  less  accuracy  than  was  available 
from  the  original  track  of  data. 

Numerous  benefits  would  result  Irom  a  more  accurate  method  of  orbit  estimation.  Also,  as  radar 
capabilities  improve  and  space  access  eases,  the  SCC  workload  will  undoubtedly  increase.  For  both  of 
these  reasons,  a  method  of  orbit  estimation  which  allows  for  some  distributed  data  processing  (at  the  sites) 
and  an  increase  in  catalog  accuracy  is  desirable. 

1.3  Method  of  Solution 

The  purpose  of  this  thesis  is  to  develop  a  means  of  compressing  the  data  in  a  track  of 
observations  while  retaining  the  majority  of  acciuacy  available  in  the  track.  Specifically,  least  squares 
differential  correction  is  used  to  determine  a  state  vector  and  associated  covariance  matrix  which  very 
efficiently  describes  all  of  the  data  in  the  original  track.  This  technique  is  applied  using  two  d^amics 
models  for  comparison.  One  model  uses  a  T^lor  series  expansion  of  the  position  vector  which  allows 
direct  calculation  of  the  position  and  state  transition  matrix  at  each  observation  time.  This  expansion 
includes  first  and  second  order  two-bodjy  acceleration  as  well  as  first  order  J2  zonal  harmonic  acceleration. 
The  second  model  uses  numerical  integration  of  the  equations  of  motion  (EOMs)  and  the  equations  of 
variation  accounting  for  two-bod^  and  J2  accelerations. 

To  demonstrate  the  accuracy  retained  from  this  procedure,  the  scope  of  the  thesis  was  extended  to 
include  orbit  estimation  using  state  vectors  and  covariance  matrices  from  multiple  tracks  of  data, 
simulating  the  combination  of  inputs  from  multiple  sites  of  the  SSN.  This  estimation  was  again 
accomplished  using  least  squares  differential  correction.  The  d^mamics  model  consisted  of  numerical 
integration  of  the  EOMs.  This  time,  the  EOMs  included  acceleration  due  to  air  drag  using  an 
atmospheric  density  model.  The  state  vector  determined  included  a  ballistic  coefiBcient,  and  is  thus 
termed  orbit  estimation.  The  associated  covariance  matrix  was  calculated  as  an  indication  of  the  accuracy 
of  the  estimate. 
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In  order  to  show  the  applicability  of  this  technique  to  the  mission  of  tracking  an  entire  catalog  of 
thousands  of  objects,  the  scope  of  the  thesis  was  again  broadened  to  include  a  one-fiftieth  scale  simulation 
of  the  tracking  and  orbit  estimation  of  the  entire  catalog.  The  justification  for  this  scale  will  be  discussed 
in  the  sixth  chapter.  Simulation.  This  simulation  included  the  satellite  scheduling,  observation,  track 
compression,  oibit  estimation  and  catalog  maintenance  required  for  the  SCC  mission  as  well  as  the 
propagation  of  the  truth  model  needed  to  generate  the  raw  data  used  as  input  observations. 

1.4  Thesis  Description 

This  thesis  is  divided  into  seven  chapters.  The  first  chapter  gives  the  background,  problem 
description  and  solution.  The  second  discusses  the  relevant  literature  and  the  theoretical  background. 

The  third  chapter  details  the  solution  methodology  and  results  of  the  Taylor  series  method  of  track 
compression  portion  of  the  thesis.  Chapter  four  details  the  solution  methodology  and  results  of  the 
integrator  method  of  track  compression.  Chapter  five  does  the  same  for  the  Global  Estimate  effort. 
Chapter  six  does  the  same  for  the  three  phases  of  the  simulation.  Finally,  chapter  seven  concludes  with 
recommendations.  To  allow  for  continuity  of  the  report,  several  detailed  derivations  are  removed  fi’om  the 
chapters  and  presented  as  appendices. 


II.  Theoretical  Background 


2.1  Introduction 

This  chapter  is  intended  to  provide  the  reader  with  the  background  in  the  area  of  batch  and 
sequential  non-linear  least  squares  orbit  determination  necessary  for  a  thorough  imderstanding  of  the 
remainder  of  the  thesis.  The  chapter  begins  with  a  brief  history  of  orbit  determination  and  progresses 
through  the  development  of  the  method  of  least  squares.  The  chapter  concludes  with  a  description  of 
modifications  to  the  method  of  least  squares  required  for  sequential  estimation. 

2  2  Orbit  Determination 

The  science  of  determining  the  orbit  of  an  object  using  observations  began  in  1609  with  Johannes 
Kepler.  He  determined  planetary  orbits  using  observations  made  the  Danish  astronomer  Tycho  Brahe. 
Kepler  determined  the  orbits  of  the  known  planets  while  formulating  his  three  laws  of  planetary  motion. 
However,  it  was  not  until  Isaac  Newton  published  The  Mathematical  Principles  of  Natural  Philosophy 
(the  Principia)  in  1687  that  the  forces  behind  Kepler’s  three  laws  were  explained.  This  began  a 
deterministic  age  in  which  the  motion  of  moons,  planets,  and  comets  could  be  predicted.  These 
predictions  were  based  on  orbits  determined  using  observation  data  consisting  of  angular  measurements  in 
the  slyf.  Boulet  and  Herget  describe  three  methods  of  orbit  determination  using  angles-only  data 
developed  by  Laplace,  Gauss  and  Olbers.  [2, 8]  Escobal  describes  an  additional  angles-only  method  and 
five  methods  using  two  position  vectors  as  well  as  several  mixed  data  techniques.  [5]  Depending  on  the 
assumptions  made  as  to  the  type  of  orbit  involved,  these  methods  require  up  to  six  independent  pieces  of 
information,  corresponding  to  the  six  constants  of  integration  obtained  fi’om  the  combination  of  Newton’s 
laws  of  acceleration  and  universal  gravitation. 

If  these  independent  pieces  of  information,  or  data,  were  perfectly  accurate,  and  the  (^mamics 
were  perfectly  known,  no  new  observations  would  be  required.  However,  in  some  cases,  there  are  forces 
acting  on  an  orbiting  bocfy,  such  as  air  drag,  which  carmot  be  exactly  predicted.  In  all  cases,  observational 
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data  contain  errors.  Thus,  there  is  always  a  need  for  more  than  the  minimum  set  of  observations.  It  is  the 
optimal  combination  of  these  imperfect  observations  which  launched  Karl  Frederich  Gauss  into  the 
history  books. 

In  1801,  Gauss  developed  the  method  of  least  squares  while  attempting  to  determine  the  orbit  of 
the  first  minor  planet  to  be  discovered,  Ceres.  [5, 13,  20]  He  later  perfected  this  method  and  published  it 
in  1809.  This  method  forms  the  basis  for  modem  estimation  theory.  It  assiunes  that  all  observations 
contain  regular  or  constant  errors  and  irregular  errors.  Gauss  states,  “it  is  up  to  the  observer  to  ferret  out 
all  sources  of  constant  error  and  remove  them.”  [6:  5]  “Inegular  errors  are  essentially  different. .  .we  have 
to  put  up  with  them  in  the  observations  themselves;  however,  we  should  reduce  their  effects  on  derived 
quantities  as  fyr  as  possible  using  judicious  combinations  of  the  observations.”  [6:  5]  The  selection  of 

this  judicious  combination  is  rooted  in  the  fundamentals  of  probability,  the  topic  of  the  next  section. 

2.3  Probability  [6,  17] 

Gauss  described  a  continuous  density  probability  function  as  a  means  of  explaining  observational 
errors.  The  probability  of  a  particular  instrument  producing  an  error  between  X]  and  X2  on  a  particular 
observation  is  obtained  by  integrating  the  instrument’s  continuous  density  probability  function,  f  (e),  over 
the  range  of  xi  to  X2. 

*2 

P(x,  <e<X2)  =  Jf(e)de  (2.1) 

*1 

A  continuous  density  probability  fimction  is  normalized. 

aa 

P(^<e<oo)=  j  f(e)de  =  l  (2.2) 

— OO 

Gauss  described  three  such  density  functions  which  will  be  repeated  here.  The  first  function  has  -a  and  +a 
as  the  limits  of  all  possible  errors,  with  all  errors  between  these  limits  being  equally  probable.  In  this 
case,  f  (e)  would  be  described  as 
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e<-a 


f(e)  = 


_L 

2a 

0 


-a<e<a 


a  <  e 


(2.3) 


This  density  function  is  shown  in  Figure  2-1. 


In  this  case,  the  mean  error,  standard  deviation  and  variance  are  given  below. 


x  =  0 


Gauss’s  second  density  function  is  described  as 


f(e)  = 


0  e<  -a 
a  +  x 


aa 

a-x 


-a<e^0 

0  <  e  <  +a 
aa 

0  a  <  e 


(2.4) 


This  density  function  is  shown  in  Figure  2-2. 
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-a  0 

Figure  2-2  Density  Distribution  Number  Two 


This  density  function  is  shown  in  Figure  2-3. 


-h  0  +h 

Figure  2-3  Density  Distribution  Number  Three 


In  this  case,  the  mean  error,  standard  deviation  and  variance  are  given  as 


24 


x  =  0 


o  =  h 


hi 

2 


Substituting  h  =  Cyfliato  equation  ( 2.5 )  returns  the  more  familiar  form  of  the  Gaussian  or  normal 


distribution. 


exp 


2  o' 


J 


(2.6) 


2  4  Central  Limit  Theorem 

Equation  ( 2.6 )  is  remarkable  for  its  frequent  occurrence  in  nature.  One  of  the  primary  reasons 
for  this  is  described  by  the  Central  Limit  Theorem.  The  Central  Limit  Theorem  is  stated  as  follows. 

If  random  samples  of  n  observations  are  drawn  from  a  population  with  finite  mean,  p, 
and  standard  deviation,  a,  then,  when  n  is  large,  the  sample  mean  will  be  approximately 
normally  distributed  with  mean  equal  to  |i  and  standard  deviation  G  /  ^^n  .  The 
approximation  will  become  more  and  more  accurate  as  n  becomes  large.  (12:  153] 

The  Central  Limit  Theorem  has  many  applications  to  the  field  of  estimation  theory.  One  of  these  is  the 

error  probability  density  function  for  an  instrument.  If  the  instrument  is  well  designed,  no  single  source  of 

error  will  dominate  the  results.  Instead,  the  irregular  errors  which  Gauss  described  will  result  from  the 

summation  of  many  small  errors,  none  of  which  dominate.  The  Central  Limit  Theorem  states  why  such 

instruments  exhibit  errors  which  follow  a  normal  distribution.  This  is  quite  significant  because  we  will 

not,  in  general,  know  the  density  functions  of  the  individual  sources  of  errors.  The  Central  Limit 

Theorem  states  that  the  resulting  combination  will  be  normal  regardless  of  the  shape  of  the  contributing 

error  density  functions.  The  only  requirement  is  that  the  individual  distributions  each  contribute  an 

infinitesimal  standard  deviation,  and  that  many  different  sources  are  combined. 

The  Central  Limit  Theorem  allowed  Gauss  to  develop  the  method  of  least  squares  by  enabling 

him  to  mathematically  describe  the  shape  of  the  error  curves  without  actually  knowing  all  of  the  sources. 

Thus,  it  is  a  cornerstone  for  this  thesis.  A  second  area  of  this  thesis  where  the  Central  Limit  Theorem  is 

used  comes  in  the  random  number  generator  used  to  add  noise  to  the  observational  data.  A  random 
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number  generator  which  produces  random  numbers  with  a  distribution  similar  to  Figure  2-1  is  called 
multiple  times,  and  the  results  are  summed.  By  the  Central  Limit  Theorem,  this  results  in  the  generation 
of  random  numbers  with  a  distribution  approximating  a  normal  distribution,  as  do  the  real-world 
observations. 

2.5  Least  Squares 

The  Central  Limit  Theorem  and  equation  ( 2.6 ),  the  normal  distribution  (also  known  as  the 
Gaussian  distribution),  are  both  fundamental  to  Gauss’s  method  of  least  squares.  To  understand  the 
development  of  Gauss’s  method,  we  must  first  describe  the  problem  that  was  at  hand.  The  minor  planet 
Ceres  had  been  discovered  and  then  lost.  There  were  more  than  enough  observations  before  it  was  lost  to 
provide  the  six  pieces  of  information  necessary  for  orbit  determination.  Astronomers  attempted  to 
reacquire  the  object  fitting  orbits  to  the  observational  data,  but  this  data  was  not  of  sufficient  accuracy. 
Gauss  theorized  that  the  best  orbit  to  be  calculated  would  not  necessary  agree  precisely  with  any  one  piece 
of  data,  but  would  instead  disagree  as  little  as  possible  with  all  of  the  data.  To  calculate  this  orbit.  Gauss 
developed  the  method  of  least  squares  in  its  multi-dimensional,  non-linear  form.  Here,  we  will  develop  it 
in  a  single  dimensional  form,  expand  it  to  its  multi-dimensional,  linear  form  and  finally  to  its  full  non¬ 
linear  form. 

2.5. 1  One-Dimensional  Least  Squares.  If  n  independent  observations  of  a  single  quantity,  x,  are 
taken  (call  this  data  Z],  Z2, ...  Zn)  with  instruments  with  standard  deviations  of  CTi,a2, ...  On,  these  will 
form  a  set  of  errors  of  ei,  e2,  ...  ^  where  ei  =  x  -  Zj.  We  can  determine  the  joint  probability  of  having 
obtained  this  particular  set  of  data  which  we  will  call  f(z),  as 


where  we  could  solve  for  this  probability  if  we  knew  the  actual  value  of  x.  But  this  value  is  unknown. 
This  is  where  Gauss  applies  the  fundamental  Principle  of  Maximum  Likelihood.  Instead  of  knowing  the 
true  value  of  x,  and  using  equation  ( 2.7 )  to  determine  the  probability  of  having  gotten  this  set  of  data. 
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Gauss  assumes  that,  because  we  have  alrea^  gotten  this  set,  the  true  value  of  x  must  be  that  value  which 
maximizes  the  probability  of  having  gotten  that  set.  Thus,  a  value  of  x  is  determined  which  drives 
equation  ( 2.7 )  to  a  maximum.  This  occurs  when  the  terms  inside  the  exponentials  are  driven  to  a 
maximum.  At  a  glance,  one  can  tell  this  equates  to  finding  the  value  of  x  which  is  closest  to  all  of  the  z; 
with  the  errors  being  squared,  thus  the  description,  the  method  of  least  squares.  The  problem  has  thus 
become  one  of  finding  a  value  for  x  (we  will  call  this  value  x )  which  minimizes  the  following  sum 


I 


(2.8) 


for  given  values  of  Zi  and  cTj.  The  value  of  x  is  obtained  by  differentiating  equation  ( 2.8 )  with  respect  to 
X  and  setting  the  result  equal  to  zero. 


z- 

~ 


=  0 


or 


X  = 


(2.9) 


The  denominator  of  equation  ( 2.9 )  provides  a  measiu^  of  the  relative  weight  of  the  estimate.  In  fact,  the 
variance  of  this  estimate  is  given  [20] 


(2.10) 


Thus,  equations  ( 2.9 )  and  (2.10)  provide  the  best  fit  estimate  of  the  true  value  of  x  and  a  statistical 
measure  of  the  reliability  of  this  estimate,  . 

2.5.2  Multi-Dimensional  Linear  Least  Squares.  [3,  7,  20]  To  extend  the  results  of  the  previous 
section  to  a  multi-dimensional  case,  the  value  to  be  estimated  must  be  redefined  as  a  state  vector,  x.  The 
individual  pieces  of  data  are  now  typically  vectors  as  well.  These  will  be  denoted  as  Zi.  These  vectors  are 
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not  necessarily  of  the  same  order  as  the  state  vector  because,  in  the  linear  least  squares  development,  we 
no  longer  assume  that  the  instnunent  directly  measures  the  values  of  the  state.  In  fact,  the  instrument  is 
not  required  to  measure  its  quantities  at  the  time  of  interest.  Instead,  a  relationship  is  developed  between 
the  state  vector  at  the  time  of  interest,  which  we  will  call  the  reference  time,  and  the  state  vector  at  any 
other  time.  This  relationship  is  based  on  the  state  transition  matrix.  The  state  transition  matrix  allows  us 
to  determine  the  value  of  the  state  vector  at  any  time  based  on  the  value  of  the  state  vector  at  the  reference 
time  and  the  elapsed  time  in  between.  Because  this  development  assumes  linearized  Ramies,  the  state 
transition  matrix,  denoted  O  (t,  ) ,  is  simply  multiplied  the  state  vector  at  the  reference  time,  x  (to  ) . 

x(t)  =  <I»(t,to)x(to)  (2.11) 

As  mentioned,  the  instrument  does  not  directly  measure  the  state  vector.  Instead,  it  measures 
some  parameters  related  to  the  state  vector  and  includes  the  ever  present  errors.  In  this  development,  this 
observation  relation  is  also  linear  and  is  expressed  as 

Zi=HiX(tj)  +  ei  (2.12) 

By  combining  the  two  equations,  we  see  the  relationship  between  the  instrument  data  and  the  state  vector 
at  the  reference  time. 

Zj=Hi<5(tj,to)x(to)  +  ei  (2.13) 

Finally,  a  covariance  matrix,  Qi,  is  now  introduced.  This  is  a  square  matrix  composed  of  the 
instrument  variances  for  each  element  of  the  observation  vector  as  the  diagonal  elements  and  the 
covariances  as  the  off-diagonals.  The  covariances  represent  the  degree  of  statistical  interdependence 
between  the  error  of  one  observation  and  the  error  of  another.  In  many  applications,  these  covariances  are 
zero,  making  the  covariance  matrix  a  diagonal  matrix. 

The  following  abbreviation  is  common. 

TiSHi4»(ti,to)  (2.14) 

Additionally,  for  N  independent  observation  vectors,  the  matrices  and  vectors  can  be  grouped  as  follows. 
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(2.15) 


z, 


[0  0  -  QnJ 

Using  the  same  method  as  section  2.5,  the  following  estimate  of  the  state  vector,  x(tQ ) ,  results.  [20] 

i(t.)=(T^Q-'TrT^Q-'Z  (2.16) 

The  reliability  of  this  estimate  is  expressed  as  the  covariance  matrix,  .  [20] 

P,  =(t^Q-'t)'’  (2.17) 

2.5.3  Non-Linear  Least  Squares.  [20]  The  non-linear  least  squares  development  differs  from  the 
linear  development  in  two  ways.  First,  the  H  matrix  and  the  O  matrix  are  obtained  by  linearizing  about  a 
reference  condition.  Second,  the  method  produces  a  correction  to  a  reference  state,  and  must  therefore  be 
iterated  until  this  correction  is  smaU. 

Instead  of  using  a  state  transition  matrix,  O  (t  j ,  ) ,  which  converts  the  state  vector  at  the 

reference  time  to  a  state  vector  at  atty  other  time,  its  linearization  is  instead  used  to  convert  a  small  change 
in  the  reference  state  to  a  small  change  in  the  state  at  another  time. 

5x(ti)  =  0(tj,to)8x(to)  (2.18) 

Unfortunately,  this  linearization  matrix  is  denoted  the  same  and  is  still  called  the  state  transition  matrix. 

A  more  appropriate  term  might  be  state  correction  transition  matrix.  The  linearization  occurs  about  a 
reference  state.  This  reference  state  is  obtained  (in  the  oibit  determination  problem)  from  an  initial  orbit 
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determination.  The  methods  described  by  Boulet,  Herget  and  Escobal  would  all  be  appropriate  for  this 
initial  reference  state. 

Instead  of  the  observation  relation 

2;  =Hi  x(ti)  +  ei  (2.12) 

the  non-linear  relation  becomes 


2;  =Gi(x(ti))+ei 

This  is  linearized  about  the  reference  state  at  the  observation  time  as  follows. 


dx 


(2.19) 


(  2.20  ) 


The  T  matrix  and  the  Q  matrix  are  defined  as  in  the  linear  case  and  a  residual  vector  is  introduced  as 

r,  =2j-Gj(x,rf(ti))  (2.21) 

The  estimate  of  the  state  vector  is  given  as  [20] 

x  =  Xref(to)  +  Sx(to)  (2.22) 

where 

5s(t<,)=(T^Q-'T)'‘T'  Q-'r 

and  the  covariance  matrix  is  [20] 

Ps  =  (t”"  Q~’ T)"'  (2.23) 

Equations  ( 2.22 )  through  ( 2.23  )  provide  the  means  of  determiiung  the  state  vector  at  the 
reference  time,  and  the  associated  covariance  matrix.  However,  there  is  one  more  important  step  in  the 
method  of  least  squares.  The  resulting  reference  state  can  be  compared  to  the  observation  data  to 
determine  the  set  of  residuals,  fi.  These  residuals  should  agree  statistically  with  the  instrument’s 
accuracies.  For  example,  if  the  instrument  is  known  to  have  a  standard  deviation  of  ten  meters  in  one 
parameter,  then  approximately  68%  of  the  residuals  should  be  less  than  ten  meters.  This  check  of  the  data 
is  required  in  order  to  show  that  something  has  not  gone  drastically  wrong  with  the  estimation  algorithm. 
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2.5.4  Computer  Implementation.  [20]  A  modification  to  the  above  equations  allows  more 
efficient  computer  implementation.  Instead  of  using  the  large,  combined  matrices  described  by  equations 
( 2. 15  ),  the  individual  Ti  and  Qi  matrices  and  the  individual  Zi  vector  can  be  used  in  the  following 
equations.  [20] 


TV  Qi  ’  T, 

V  i 


(2.24) 


These  equations  are  iterated  until  5  x  has  converged  to  a  small  enough  value.  The  criteria  for 
convergence  can  come  from  the  covariance  matrix.  If  the  correction  to  each  element  of  the  state  vector  is 
significantly  smaller  than  the  square  root  of  the  corresponding  diagonal  element  of  the  covariance  matrix, 
there  is  no  need  to  continue  iterating. 


2. 6  Sequential  Estimation 

The  methods  described  in  section  2.5  are  termed  batch  estimators.  They  take  all  of  the  available 
data  and  process  it  at  one  time.  This  is  the  most  accurate  means  of  data  reduction  but  can  become 
cumbersome  for  certain  i^ems.  In  the  continual  process  of  orbital  estimation  for  a  near-Earth  satellite, 
data  may  be  available  spanning  several  decades.  It  would  be  too  time  intensive  to  attempt  to  process  all  of 
this  data  to  produce  the  latest  orbit  estimate.  If  the  epoch  time  of  the  estimate  is  to  be  kept  somewhat 
ciuTent,  the  (fynamics  must  be  integrated  from  the  time  of  the  older  data  to  the  current  epoch.  This  can 
take  considerable  time  and  is  not  entirely  accurate  since  the  <^iiamics  model  used  is  never  perfectly 
accurate.  Instead,  a  common  practice  is  to  estimate  an  orbit  using  the  new  data  since  the  last  orbit 
estimate  along  with  the  results  of  the  previous  estimate.  Thus,  the  previous  estimate  acts  like  data  to  the 
new  estimate  and  the  previous  covariance  matrix  acts  like  an  observation  covariance  matrix.  This  process 
is  repeated  as  often  as  necessary.  Thus,  each  estimate  is  at  a  current  epoch  and  the  older  data  is  included 
in  the  process  by  including  the  previous  estimate  (which  was  based  on  the  older  data).  This  process  of 
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continually  discarding  old  data  and  keeping  the  estimate  and  covariance  matrix  is  called  sequential 
estimation.  The  process  of  sequential  estimation  using  the  method  of  least  squares  for  each  estimate  is 
called  a  Bayes  Filter. 

One  common  problem  with  any  sequential  estimator  is  the  reducing  size  of  the  covariance 
matrix.  Although  this  seems  like  a  desirable  result,  it  can  eventually  turn  bad.  As  the  amoimt  of  data 
going  into  the  estimate  increases,  the  covariance  matrix  will  eventually  become  so  small  that  the  estimator 
essentially  believes  it  has  achieved  perfection.  New  data  is  evaluated  but  its  observation  covariance  matrix 
will  be  quite  large  compared  to  the  covariance  of  the  previous  estimate  and  so  the  new  data  is  almost 
completely  ignored.  If  the  (fynamics  were  perfectly  known,  and  there  were  no  numerical  enors  in  the 
implementation,  this  would  be  acceptable.  However,  as  mentioned  previously,  there  are  actually  some 
influences  on  the  motion  of  the  satellite  which  caimot  be  entirely  predicted  Effects  such  as  air  drag  are 
not  entirely  predictable  because  the  atmospheric  density  at  any  given  point  is  highly  unstable.  Predictions 
may  be  only  accurate  to  an  order  of  magnitude  in  some  cases.  Thus,  the  newer  data  is  necessary  to 
provide  a  current  orbit  estimate,  and  some  means  of  preventing  the  covariance  matrix  from  shrinking 
must  be  implemented  The  method  used  in  this  thesis  is  called  fading  memory.  Fading  memory  involves 
multiplying  the  covariance  matrix  by  another  matrix,  called  a  ^-matrix,  whose  size  is  based  on  the  age  of 
the  covariance  matrix.  The  p-matrix  is  a  diagonal  matrix  whose  elements  are  made  from  the  reciprocal  of 
coefficients  raised  to  the  power  of  the  number  of  days  between  the  epoch  of  the  covariance  matrix  and  the 
epoch  of  the  current  estimate.  The  coefficients  are  typically  less  than  one  and  represent  the  time  scale  for 
the  accuracy  of  that  particular  state  vector  element. 

P„ew=P..dP’^‘  (2.25) 

Thus,  if  the  p  coefficients  are  0.933,  and  the  covariance  matrix  is  10  days  old  the  covariance  elements 
would  be  doubled  corresponding  to  a  data  half-life  of  ten  days.  Different  state  vector  elements  can  have 
different  p  coefficients  corresponding  to  different  time  scales  for  the  length  of  validity  of  the  covariance 
elements. 
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2.7  Summary 

This  chapter  has  described  the  basics  of  orbit  determination,  least  squares  batch  estimation,  the 
Bayes  Filter  and  fading  memoiy.  The  thesis  makes  extensive  use  of  the  non-linear,  multi-dimensional 
method  of  least  squares.  It  uses  the  computer  implementation  method  described  in  section  2.5.4.  The 
final  and  most  important  portion  of  the  thesis  uses  a  Bayes  Filter  with  fading  memoiy  in  a  long-term 
simulation.  The  next  chapter  begins  the  detailed  description  of  the  thesis  effort. 
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///.  Taylor  Series  Track  Compression 


3. 1  Introduction 

This  chapter  describes  the  Taylor  series  method  of  track  compression.  It  first  details  the 
development  of  the  truth  model  used  to  generate  the  various  sets  of  input  observation  data.  It  then  details 
how  the  Taylor  series  method  of  track  compression  was  developed  using  the  method  of  least  squares. 
Finally,  it  presents  the  results  of  the  Taylor  series  track  compression  portion  of  the  thesis. 

3.1.1  Notation.  Throughout  the  thesis,  bold-faced  letters  will  be  used  to  represent  vectors.  State 
vectors  are  represented  as  X .  Position  and  velocity  vectors  are  represented  as  R  and  V .  Components  of 
position  and  velocity  will  be  represented  as 


Table  3-1  Notation 


position 

velocity 

[Ri 

R, 

RkF 

[V, 

V, 

v.r 

K 

Ry 

[Vx 

Vv 

v.r 

[R. 

R, 

[V, 

V2 

vF 

y 

3. 1.2  Units.  To  a  large  extent,  canonical  units  have  been  used  throughout  the  thesis.  Distances 
are  measured  in  Distance  Units  (DUs),  time  is  measured  in  Time  Units  (TUs),  seconds  or  Julian  Days 
(JD).  A  DU  is  defined  as  the  mean  equatorial  radius  of  the  Earth  equal  to  6378. 137  km  [22].  A  TU  is 
defined  as  the  time  that  a  satellite  in  a  circular,  Keplerian  orbit  with  a  semi-major  axis  of  one  DU  would 
take  to  travel  a  distance  of  one  DU.  This  is  equal  to  13.446851 15881  minutes  [22].  Exceptions  to  the 
canonical  units  standard  occin  in  the  time  steps  used  the  numerical  integrator  and  the  reporting  of 
observation  times. 
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3.2  Truth  Model 


The  truth  data  is  generated  using  a  foiulh-order  Runge-Kutta  numerical  integrator.  This 
integrates  a  set  of  first-order  differential  equations  of  the  form 

dx 

-  =  f(x)  ,3.1, 

The  formulas  for  the  standard  fourth-order  Runge-Kutta  integrator  are 

^n+i  =  Xjj  +  ~  (kj  +  2  kj  +  2  k3  +  k4  )  ( 3.2 ) 

where  kj  =  At  f (x^) 

kj  =Atf^x. +^j 
=  At  f(x„ -I- kj)  [1, 11] 

In  this  application,  the  variable  of  interest  is  a  state  vector,  X ,  and  the  differential  equation  takes  the 
vector  form 

dX 

—  =  f(X)  (3.3) 

3.2. 1  State  Vector.  For  the  track  compression  portions  of  the  thesis,  the  state  vector,  X ,  is 
defined  as 

x=[r,  R,  R,,  V,  V,  (3,4) 

3.2.2  Equations  of  Motion.  For  the  track  compression  portions,  the  EOMs  are  derived  from  the 
Earth’s  geopotential  accounting  for  the  J2  zonal  harmonic  caused  ly  the  Earth’s  oblateness.  Escobal  gives 
this  potential,  <1> ,  (not  to  be  confused  with  the  state  transition  matrix  with  the  same  symbol),  as  [5] 
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where  G  is  the  gravitational  constant 
m  is  the  mass  of  the  Earth 

is  the  equatorial  radius  of  the  Earth 
r  is  the  distance  from  the  Earth’s  center 
h  =  1082.28  ±  0.3  X  10-®  (unitless) 
sin  8  =  z/t 

The  acceleration  due  to  gravity  about  a  planet  is  given  as  the  gradient  of  the  potential,  or 

-^  =  VO  (3.6) 

and  these  terms  are  given  as 

f  =  -^[l+f^(l-5sin-S)]  (3.7) 

dy  3x  J 

BO  Gmz  r  •,.,2?) 

-^  =  — ^|^l  +  --|j-(3-5sin  5j 

Using  canonical  units  aHows  the  substitution  |i  =  Gm  =  1  and  R®  =  1.  The  equations  of  motion  become 
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where  r  is  detenmned  from 


r  =  Vx'+y'+z^  (3.9) 

3.2.3  Reference  Frames.  The  truth  model  tracks  the  state  vector  and  acceleration  components  in 
the  Earth-C^entered,  Inertial  (ECI)  reference  frame.  This  is  sometimes  referred  to  as  the  IJK  frame,  as 
these  letters  represent  the  three  basis  vectors.  This  frame  uses  the  center  of  the  Earth  as  its  origin,  the 
equatorial  plane  as  the  fundamental  frame,  the  vernal  equinox  as  the  direction  of  the  first  basis  vector,  and 
the  North  pole  as  the  direction  of  the  third  vector.  The  second  vector  is  in  (he  equatorial  plane, 
perpendicular  to  the  first  and  third  vectors  in  a  right-handed  sense. 

An  Earth-Centered,  Rotating  (ECR)  fiame  is  also  used.  This  frame  is  similar  to  the  ECI  frame 
except  that  the  first  basis  vector  always  points  out  from  the  center  of  the  Earth  through  the  equator  at  a 
longitude  of  zero  degrees  (the  longitude  of  Greenwich).  To  speed  generation  of  observation  data,  the  site 
position  vector  is  determined  once  in  the  ECR  fiame.  The  state  vector  of  the  satellite  is  then  rotated  to  the 
ECR  fiame  at  each  time  step  and  compared  to  the  site  vector  to  determine  visibility.  By  rotating  the  state 
vector  through  a  single,  simple  rotation  about  the  third  axis,  the  site  vector  does  not  need  to  be 
recalculated  at  each  time  step.  It  remains  constant  in  the  ECR  fiame. 

Finally,  the  Topocentric-Horizon  reference  fiame  is  used  to  determine  the  range,  azimuth  and 
elevation  of  the  satellite  as  seen  fiom  the  site.  This  reference  fiame  has  the  site  as  its  origin,  the  local 
horizon  as  the  fundamental  plane.  South  (in  the  local  horizon)  as  the  first  direction.  East  as  the  second 
direction,  and  straight  up  (the  zenith)  as  the  third.  This  fiame  is  usually  referred  to  as  the  SEZ  fiame 
owing  to  the  direction  of  the  basis  vectors.  The  SEZ  frame  allows  easy  calculation  of  the  range,  azimuth 
and  elevation  observations  which  would  be  reported  Ity  the  sites. 

3.2.4  Observation  Data.  The  truth  model  is  used  to  generate  the  observation  data  for  use  by  the 
Taylor  series  track  compression  algorithm.  The  observation  data  represents  a  typical  set  of  observations 
which  a  site  of  the  SSN  could  generate  as  a  satellite  passes  within  its  field  of  view.  The  radar  sites  used 
will  be  described  next. 
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3.2.4. 1  Radar  Locations.  The  radar  site  locations  shown  in  Table  3-2  were  used  for  this 
thesis.  These  are  the  actual  SSN  locations.  [19] 

Table  3-2  Radar  Sites  of  the  SSN 

Station  Latitude  (deg)  Longitude  (deg)  Altitude  (m) 
liiS  -4.67174786  55.47782059  560.50 

Reef  -7.27003056  72.36999860  -68.375 

Guam  13.61518782  144.85604938  218.93 

Hula  21.56226524  201.75789406  429.42 

Cook  34.82259890  239.49814705  271.53 

38.805943055  255.471532222  1899.42 

loss  42.94782144  288.37343743  203.37 

Pogo  76.51536439  291.40114169  147.03 

Lion  51.11758338  359.0936545  146.59 

The  world-wide  distribution  of  these  sites  is  evident  fiom  Figure  3-1  below. 


Figure  3-1  Radar  Sites  of  the  SSN 


3. 2.4. 2  Radar  Capabilities.  Observation  data  is  generated  assuming  that  all  sites  have 
the  same  type  of  radar.  These  radar  measure  the  range,  azimuth  and  elevation  of  the  satellite.  The  rate  of 
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observations  varied  for  different  portions  of  the  thesis.  For  the  track  compression  development,  data  rates 
varied  from  one  observation  every  five  seconds  to  as  high  as  twelve  observations  per  second. 

All  radar  sites  were  modeled  with  the  same  accuracy.  The  standard  deviations  used  were  100 
meters  in  range,  0.025  <fegrees  in  azimuth  and  0.025  degrees  in  elevation.  These  are  typical  SSN 
accuracies.  [19] 

3. 2. 4. 3  ECR  Site  Vector  Calculation.  To  calculate  the  position  vector  of  each  site  in  the 
ECR  frame,  one  first  needs  to  understand  the  effects  of  the  Earth’s  oblateness  on  latitude,  and  the 
distinction  between  geocentric  and  geodetic  latitudes.  As  shown  in  Figure  3-2,  geocentric  latitude,  L',  is 
the  angle  between  the  equatorial  plane  and  a  line  from  the  site  to  the  center  of  the  Earth.  Geodetic 
latitude,  L,  is  the  angle  between  the  equatorial  plane  and  a  line  through  the  site  and  perpendicular  to  the 
local  horizon.  This  is  what  is  meant  when  latitude  is  typically  reported. 


The  latitudes  in  Table  3-2  are  geodetic  latitudes.  To  determine  the  ECR  position  vector  of  a  site  with 
geodetic  latitude  L,  longitude  X,  and  altitude  H,  the  following  equations  are  used.  [1] 


R 


X  cosX, 
X  sin  A. 
z 


(3.10) 
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Vl-e^  sin^  L 
ae(l-e') 


Vl-e^sin^ 


+  H 


+  H 


cosL 


sinL 


where  Be  is  the  mean  equatorial  radius  of  the  Earth  (1  DU) 

e^  is  the  squared  eccentricity  of  the  Earth  (the  amount  of  oblateness) 

From  [22]  e^  =  0.00669437999013. 

3. 2.4.4  Observation  Geometry.  To  predict  the  measurements,  the  position  vector  of  the 
satellite  relative  to  the  site  is  first  calculated  in  ECR  coordinates  and  then  rotated  into  the  SEZ  frame. 

The  position  components  of  the  state  vector  are  used  to  create  the  vector  of  the  position  of  the  satellite, 
relative  to  the  center  of  the  Earth,  expressed  in  ECl  coordinates.  This  vector  is  then  rotated  about  the 
third  axis  an  angle  equal  to  the  local  sirtereal  time  at  Greenwich,  resulting  in  a  position  vector  in  ECR 
coordinates.  The  position  vector  of  the  site  is  then  subtracted  from  this  to  leave  the  position  vector  of  the 
satellite,  relative  to  the  site,  expressed  in  ECR  coordinates.  This  vector  is  then  rotated  in  a  positive 
direction  about  the  third  axis  1^  the  longitude  of  the  site,  and  then  in  a  negative  direction  about  the  second 
axis  the  site  colatitude,  (equal  to  90°  minus  the  latitude).  This  results  in  a  vector  representing  the 
position  of  the  satellite,  relative  to  the  site,  expressed  in  the  topocentric-horizon  reference  frame.  This 
vector  is  called  RhoSEZ  and  is  abbreviated  - 

Figure  3-3  Observation  Geometry,  on  the  next  page,  shows  the  topocentric-horizon  reference 
frame  and  the  geometry  involved  in  determining  the  radar  observations. 
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Figure  3-3  Observation  Geometry 

This  figure  shows  the  position  vector  of  the  satellite  relative  to  the  site.  This  vector  is  converted  to 
measurements  of  range,  azimuth  and  elevation  using  the  following  equations. 

rstigc  —  IftEZ  I  ~  a/Ps  Pe  Pz  (3.11) 

Pe 

azimuth  =  -  tan  — 

Ps 

Pz 

elevation  =  tan  ^  —==== 

Vps'+Pe' 

3.2.5  Noise.  Equations  (3.11)  provide  range,  azimuth  and  elevation  data  for  processing  by  the 
track  compression  algorithms.  This  is  not,  however,  realistic  data.  As  discussed  in  chapter  two,  real  data 
consists  of  the  true  portion,  and  the  always  present,  never  separable  error,  the  noise.  Observations  of  the 
SSN  contain  contributions  fi'om  many,  very  small,  error  sources.  These  errors  combine  to  form  the 
reported  accuracy  of  the  instrument.  According  to  the  Central  Limit  Theorem  (see  chapter  two),  the 
errors  for  these  measurements  follow  a  Gaussian  distribution.  This  means  that  the  perfect  data  obtained 
fi'om  equations  (3.11)  must  be  intentionally  altered  to  obtain  noisy  data.  This  process  occurs  by  adding 
the  product  of  a  random  number  and  the  standard  deviation  for  that  ^pe  of  data.  The  random  number  is 
provided  by  a  routine  which  determines  random  numbers  of  a  Gaussian  distribution  with  a  mean  value  of 
zero  and  a  variance  of  one.  Figure  3-4  shows  the  distribution  of  10,000  of  these  random  numbers. 
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standard  deviations 


Figure  3-4  Gaussian  Random  Number  Distribution 


This  error  closely  resembles  a  Gauss  curve  with  a  variance  equal  to  one.  The  following  three  figures  show 
the  ou^t  of  noisy  data  from  a  track  of  100  observations.  Horizontal  lines  indicate  standard  deviation 
intervals  on  each  figure.  These  show  the  effect  of  adding  noise  to  the  observations  for  realism. 


Figure  3-5  Random  Range  Noise 
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elevation  noise  (degrees)  azimuth  noise  (degrees) 


0.075 


4  ♦ 


Figure  3-6  Random  Azimuth  Noise 


Figure  3-7  Random  Elevation  Noise 

3. 3  General  Description  of  Track  Compressi  on 

The  two  methods  of  track  compression  developed  in  this  thesis  are  based  on  the  use  of  the  multi¬ 
dimensional,  non-linear  least  squares  algorithm,  described  in  chapter  two,  to  determine  a  state  vector 
which  represents  the  actual  state  of  the  satellite  at  a  particular  time.  The  method  also  produces  the  very 
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important  covariance  matrix.  The  algorithm  essentially  says  that  if  a  satellite  were  at  the  state  computed 
at  the  reference  time,  then  it  would  describe  an  arc  across  the  sky  very  close  to  the  arc  actually  observed  by 
the  site.  In  addition,  the  covariance  matrix  describes  the  accuracy  of  the  estimated  state  vector. 

The  state  computed  at  each  iteration  of  the  least  squares  algorithm  is  called  the  reference  state. 
The  time  of  this  state  is  the  reference  time.  The  trajectory  which  the  satellite  would  trace  if  it  were 
actually  at  the  reference  state  at  the  reference  time  is  called  the  reference  trajectory.  The  reference 
trajectory  is  used  to  predict  observations  to  compare  with  the  actual  observations  to  determine  the 
accuracy  of  the  reference  state. 

This  method  does  not  quite  determine  the  satellite’s  orbit,  because  it  does  not  provide 
information  about  the  important  drag  parameter,  the  ballistic  coefiBcient.  The  time  scale  required  for  air 
drag  to  take  effect  is  too  large  to  be  seen  during  the  relatively  short  arcs  observed  as  one  track  of  data. 
Also,  the  information  is  only  based  on  observations  taken  on  a  small  piece  of  the  orbit  (it  excludes  input 
from  the  opposite  side  of  the  orbit).  If  the  data  were  perfectly  accurate,  this  arc  could  be  extrapolated  to 
completely  describe  the  entire  orbit.  However,  since  the  information  determined  by  the  method  is  only  an 
estimation,  the  errors  involved  in  predicting  the  satellite’s  position  during  portions  of  the  orbit  far  from 
the  arc  will  be  greatly  increased  For  these  reasons,  while  the  process  of  track  compression  does  produce  a 
six-element  state  vector  which  satisfies  the  classic  definition  of  orbit  determination  [5],  for  purposes  of 
this  thesis,  the  actual  orbit  estimation  is  not  considered  to  occur  until  results  from  multiple  tracks  are 
combined  (see  chapter  five). 

The  two  methods  both  require  a  dynamics  model  to  predict  what  observations  would  be  seen  if 
the  reference  state  were  the  actual  state.  This  (fynamics  model  is  used  to  take  a  reference  state  at  the 
reference  time  and  compute  the  state  of  the  satellite  at  each  observation  time.  These  states  are  then  used 
to  predict  the  observations.  The  two  methods  of  track  compression  developed  in  the  thesis  differ  in  the 
type  of  dynamics  model  used. 
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3.4  Detailed  Taylor  Series  Track  Compression  Development 

The  first  method  of  track  compression  developed  for  this  thesis  involves  a  simplified  dynamics 
model  using  a  Taylor  series  expansion  to  determine  the  position  vector  of  a  reference  orbit  at  each 
observation  time.  This  description  disagrees  slightly  with  that  of  the  preceding  section.  Typically,  in  the 
method  of  non-linear  least  squares,  the  entire  state  is  ^termined  at  each  observation  time.  In  this 
application,  however,  only  the  position  need  be  computed.  The  observation  data  do  not  include  any  rate 
information.  Thus,  the  (fynamics  model  does  not  need  to  provide  the  entire  state  at  each  time,  just  the 
position. 

3. 4. 1  Initial  Reference  State.  To  implement  the  non-linear  least  squares  algorithm  (see  chapter 
two),  a  method  of  determining  an  initial  reference  state  vector  must  first  be  developed.  This  state  must  be 
accurate  enough  to  allow  the  least  squares  algorithm  to  converge.  Additionally,  the  method  must  be  fast 
to  allow  for  computer-execution  thousands  of  times  per  day  at  each  site.  In  this  application,  the  reference 
state  is  determined  from  observation  data  near  the  middle  of  the  arc.  Specifically,  the  reference  time  is 
determined  as  the  time  midway  between  the  first  and  last  observation  in  the  arc.  Two  observations,  one 
approximately  five  seconds  before  the  reference  time  and  one  approximately  five  seconds  after  the 
reference  time  are  extracted  fi'om  the  data  set.  Each  observation  is  used  to  create  a  position  vector  at  that 
observation  time.  This  is  created  in  a  process  exactly  reversed  from  that  described  in  section  3. 2.4.4. 

A  short,  second-order  arc  is  then  fit  between  the  two  vectors;  the  reference  time  (very  near  the 
middle  of  the  short  arc)  is  used  to  interpolate  along  the  arc,  between  the  two  position  vectors,  to  determine 
the  position  vector  at  the  reference  time.  The  velocity  vector  is  determined  by  assuming  a  straight  line 
path  is  traveled  between  the  two  position  vectors  in  the  (roughly)  ten  seconds  between  the  two  observation 
times.  The  position  and  velocity  vectors  are  combined  to  create  a  reference  state  vector. 

More  accurate  methods  of  preliminary  orbit  determination  were  explored.  For  instance,  Escobal 
details  five  different  methods  of  orbit  determination  from  two  position  vectors  and  time.  In  particular,  the 
Gaussian  Iteration  method  is  well  suited  to  a  time  increment  of  roughly  ten  seconds.  Each  of  these 
methods,  while  more  accurate  than  the  arc  and  line  method  described  above,  take  considerably  longer  to 
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implement.  Due  to  the  inaccuracies  in  the  observations,  the  velocity  computed  via  the  arc  and  line 
method  could  contain  significant  error.  Nevertheless,  the  least  squares  iteration  algorithm  has  been 
shown  to  correct  errors  in  velocity  as  high  as  two  orders  of  magnitude  in  only  two  iterations.  Therefore, 
even  in  the  worst  cases,  the  arc  and  line  method  allows  convergence  faster  than  would  be  possible  using  a 
more  exact  method.  It  allows  rapid  reference  state  calculation  from  the  data  arc  with  sufficient  accuracy 
for  least  squares  convergence  in  well  over  99.9%  of  the  arcs  generated. 

3.4.2  Dynamics  Model.  After  determining  an  initial  reference  state,  a  dynamics  model  must  be 
used  to  determine  the  reference  trajectory  and  the  state  transition  matrix  at  each  observation  time.  In  the 
Taylor  series  method  of  track  compression,  a  Taylor  series  expansion  of  the  position  vector  at  the 
reference  time  is  used  to  allow  direct  calculation  of  the  position  vector  at  any  other  time  (as  opposed  to 
numerical  integration).  The  state  transition  matrix  at  an  observation  time  is  also  directly  calculated  by 
taking  the  partial  derivatives  of  each  component  of  the  position  vector  at  the  observation  time  with  respect 
to  each  component  of  the  reference  state. 

The  Taylor  series  representation  of  the  position  from  the  reference  state  is  shown  below  as 
equation  ( 3.12 ). 

R(t)  =  R(to)  +  R(to)At  +  |R(t„)At'+^R(to)At^+.-  (3.12) 

where  R(to  )  is  the  position  at  the  reference  time 
R(to)  is  the  velocity  at  the  reference  time 
R(tQ  )  is  the  acceleration  at  the  reference  time 
R(to)  is  the  time  rate  of  change  of  the  acceleration 
A  t  is  the  time  between  the  observation  time,  t,  and  the  reference  time,  to 
The  position  and  velocity  vectors  are  contained  in  the  reference  state.  The  acceleration  and  the 
acceleration  derivatives  are  obtained  from  the  equations  of  motion  for  a  satellite. 
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3.4,2, 1  Order  of  Terms.  Accelerations  resulting  from  different  forces  (i.e.  different 
zonal  harmonics)  are  calculated  separately  and  their  results  are  summed.  This  allows  the  determination  of 
the  relative  contribution  of  each  term.  Analysis  of  the  contributions  is  used  to  determine  how  many  terms 
are  required  in  order  to  maintain  a  desired  accuracy  in  the  position  vector  calculation.  Because  the 
accuracy  of  the  expansion  decreases  as  the  time  interval  increases,  the  track  of  data  must  be  split  up  into 
smaller  arcs  of  data  and  each  arc  compressed  individually.  As  more  acceleration  terms  are  included  in  the 
expansion,  the  resulting  position  vectors  are  more  accurate,  which  allows  a  track  to  be  split  into  fewer 
arcs,  but  evaluation  of  the  equations  takes  longer.  As  the  accuracy  is  decreased,  the  simpler  equations  are 
evaluated  faster,  but  the  track  must  be  split  up  into  more  arcs,  thus  requiring  the  compression  scheme  to 
be  applied  to  more  arcs.  Thus  a  tradeoff  exists  between  including  more  terms  in  the  model  and  splitting 
the  tracks  into  more  arcs. 

To  determine  the  number  of  terms  required  and  the  size  of  each  arc,  the  contribution  of  each  of 
six  terms  was  graphed  as  a  function  of  the  time  interval.  These  six  terms  are  the  first,  second  and  third 
order  two-body  terms,  first  and  second  order  h  terms,  and  the  first  order  J3  term.  This  is  included  as 
Figure  3-8. 


Figure  3-8  Contribution  of  Taylor  Series  Terms 
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Based  on  this  figure,  first  and  second  order  two-bocfy  accelerations  and  first  order  J2  accelerations  are 
included  in  the  equations.  This  allows  for  arc  lengths  of  90  seconds  (45  seconds  on  either  side  of  the 
reference  time)  while  retaining  every  term  which  could  contribute  more  than  about  three  meters  to  the 
position  of  the  satellite.  Satellites  and  debris  can  be  more  than  three  meters  in  size,  so  including  more 
terms  wotild  necessitate  accounting  for  which  portion  of  the  satellite  bo^  reflected  the  majority  of  the 
radar. 


3.4.2.2  Reference  Trajectory.  Using  equation  ( 3.12 )  and  the  reference  state,  the 
position  of  the  satellite  can  be  directly  calculated,  with  three  meter  accuracy,  at  any  observation  time 
within  45  seconds  of  the  reference  time.  This  direct  calculation  avoids  the  requirement  to  numerically 
integrate  the  equations  of  motion  between  observation  times.  The  Taylor  series  method  was  originally 
chosen  because  it  was  believed  that  this  would  save  processing  time  by  eliminating  the  requirement  to 
numerically  integrate  a  state  vector  through  many  small  time  steps  to  proceed  from  one  observation  to  the 
next.  However,  as  will  be  discussed  in  chapter  four,  this  time  savings  did  not  occur. 

The  following  is  the  complete  development  of  the  terms  of  the  components  of  equation  ( 3. 12 )  in 
terms  of  the  reference  state,  X .  To  simplify  the  presentation  of  the  equations,  the  following  notation  will 
beused. 


X  = 


R 


=  [r,  R,  R3  V.  V,  Y,f 


Ri' 

R(t.)  = 

R2 

ii 

0 

> 

II 

0 

V2 

Rs. 

V3. 

(3.13) 


(3.14) 
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R,  1-5-^ 


ROo)"  J.3  R2  2r* 


f  R,  1-5- 


R3  3-5-f- 
^  r 


3R,--V. 

R(to)  =  :|^  3R,^-V, 


3R3  —  V, 

3  r 


where  r  =  -y/Rj^  +  R2^  +  R3 


_R,V,+R2V2+R3V3 

dt*^  r 


Combining  equations  (3.12)  and  (3.14)  yields  the  following  final  equation  for  the  position  of  the 
satellite  as  a  function  of  the  reference  state  and  the  time  interval  between  the  reference  time  and  the  time 
of  the  observation,  At. 


^  R,  ^  ,  SJjRj,  R3  1*  ,  1  f3R,f 

R,  +  V,  A,-5^At’  Ja,=  V, jAt 

R=  (3.15) 

^  ^  2P  4r’  I  rM  6rH  r  J 


R3  2  3J,R3f  R3  1,  2  1  (3R3f  „'\^3 

R,  +  V,  At-^At'  -  3-5-^  At^  +^\  — ^ — V,  |At' 


6rH  r 


This  equation  applies  equally  for  time  before  the  reference  time  (where  At  is  negative)  as  well  as  after. 
The  actual  computer  code  which  implements  this  equation  takes  advantage  of  the  fact  that  each  term 
changes  only  in  the  time  interval  when  evaluating  the  position  at  dijfferent  observation  times.  To  speed 


processing,  each  term  is  completely  evaluated  with  all  but  the  At  factor  at  the  beginning  of  each  iteration 
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of  the  least  squares  algorithm.  The  position  vector  at  each  time  is  then  determined  by  multiplying  each 
term  the  appropriate  At  factor  and  summing.  The  process  of  evaluating  the  reference  trajectoiy 
consists  of  determining  the  position  vector  at  each  observation  time. 

3.4.2.3  State  Transition  Matrix.  In  addition  to  a  method  of  computing  the  reference 
trajectoiy,  the  non-linear  least  squares  algorithm  also  requires  a  method  of  computing  the  state  transition 
matrix,  ^ .  This  is  derived  from  the  dynamics  model  as  well.  In  the  Taylor  series  method,  the  state 
transition  matrix  is  developed  ly  taking  the  partial  derivatives  of  the  position  vector  (from  the  Taylor 
series  expansion)  with  respect  to  the  reference  state  vector. 

‘aR,(t)  aR,(t)  aR,(t)  aR,(t)  aR,(t)  aR,(t)' 

ax,(o)  ax^co)  3X3(0)  ax^co)  ax,(o)  ax,(o) 

aR(t)_  aR^ct)  aR^ct)  aR3(t)  dR,(t)  aR3(t) 

ax(o)  ax,(o)  3X3(0)  3X3(0)  3X4(0)  axj(o)  3x^(0)  ^  ^ 

aR3(t)  aR3(t)  aR3(t)  aR3(t)  aR3(t)  aR3(t) 

_ax,(0)  3X3(0)  3X3(0)  3X4(0)  3X,(0)  3X^(0) 

where  R(t)  is  given  ly  equation  (3.15)  and  X(0)  is  the  state  vector  which  is  related  to  the  terms  of 
equation  (3.16)  through  equation  (3.13).  The  r  and  r  terms  are  each  functions  of  the  state,  thus  the 
development  of  this  matrix  requires  multiple  uses  of  the  chain  rule  and  becomes  quite  lengthy.  This 
development  is  presented  in  its  entirety  as  Appendix  C. 

3.4.3  Observation  Relation.  The  next  step  in  the  least  squares  algorithm  implementation  is  to 
determine  the  observation  relation.  The  observation  relation,  denoted  G ,  is  usually  a  function  of  the  state 
vector  at  the  observation  time.  In  this  case,  however,  the  relation  is  only  a  function  of  the  position 
components  of  the  state.  Velocity  information  is  neither  provided  by  the  dynamics  model  nor  needed  by 
the  observation  relation.  Therefore,  the  set  of  equations  comprising  the  observation  relation  will  only  be 
functions  of  the  position,  R(t) ,  and  of  time. 

range(R,t) 

G=  azimuth(R,t)  (3.17) 

elevation(R,t) 
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The  three  quantities  observed  are  range,  azimuth  and  elevation.  These  three  quantities  are 
determined  by  a  method  similar  to  that  described  in  section  3.2.4.4  ending  with  equation  (3.11).  The 
process  is  detailed  in  the  following  three  sections. 

3.4.3. 1  Determining  the  RhoIJK  Vector.  First,  the  position  vector  of  the  site  is 
subtracted  from  the  position  vector  of  the  satellite,  with  both  vectors  expressed  in  the  ECI  coordinate 
frame.  The  resulting  vector  is  called  the  RhoIJK  vector,  also  denoted  Pp^ .  The  satellite  position  vector 
is  obtained  in  the  ECI  frame  from  equation  (3.15).  The  site  position  vector  in  the  ECI  frame,  ,  is 
determined  from  the  site’s  latitude,  L,  longitude,  X,  and  altitude,  H,  and  from  the  observation  time 
(expressed  as  a  Julian  Date)  using  equations  very  similar  to  equations  ( 3. 10 ).  [1] 

PbK=R(‘)-»d.  <3.18) 


R 


site 


X  COS0 
X  sin  6 
z 


(3.19) 


X  = 


z  = 


Vl-e^  sin^  L 


+  H 


Vl-e^  sin^  L 


cosL 


sinL 


e  =  GST+A, 

where,  as  in  equation  ( 3. 10 ),  a*  is  the  mean  equatorial  radius  of  the  Earth  and  e^  is  the  squared 
eccentricity  of  the  Earth.  The  Greenwich  Sidereal  Time,  GST,  is  determined  from  the  Julian  Date  using 
the  following  equations.  [18] 
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GST  =  1.75368559  +  628.3319705  JC 

+  6.770708127x10"*  JC"'  (3.20) 

+  co^  (JD  -  Int(ID)  -  0.5) 


where  (0®  is  the  rotational  rate  of  the  Earth,  (6.30038809866574  radians  per  day)  [22] 


JD  is  the  Julian  Date  (not  modified) 

Int  (JD)  is  the  integer  portion  of  the  Julian  Date 
JC  is  the  Julian  (Century  given  by  the  following  [18] 


JC  = 


Int(JD)  + 03 -2451545 
36525 


(3.21) 


Equation  (3.19)  gives  the  ECI  position  vector  of  the  site.  This  is  subtracted  from  the  satellite’s  ECl 
position  vector,  given  by  equation  (3.15),  yielding  the  position  vector  of  the  satellite  relative  to  the  site  in 
the  ECI  frame,  which  is  RhoIJK. 

3. 4. 3. 2  Determining  the  RhoSEZ  Vector.  The  next  step  in  determining  the  observation 
relation,  G ,  is  to  calculate  the  position  vector  of  the  satellite  relative  to  the  site  measured  in  the 
Topocentric-Horizon  (SEZ)  reference  frame.  This  reference  flume  conversion  simply  requires  two 
rotations  of  the  RhoIJK  vector.  This  is  accomplished  via  the  following  rotation  matrix. 


cos(LST)  cos(colat) 
-sin(LST) 
cos(LST)  sin(colat) 


sin(LST)  cos(colat)  -  sin(colat) 
cos(LST)  0 

sin(LST)  sin(colat)  cos(colat) 


(  3.22  ) 


where  colat  is  the  colatitude  of  the  site,  equal  to  90°  -  L.  Thus,  the  RhoSEZ  vector  is  computed  from 

ftEZ=DPDK  (3-23) 

3.4.3.3  Final  Observation  Relation  Determination.  The  range,  azimuth  and  elevation 
are  then  computed  from 


range  —  +  Pe  Pz 


(3.11) 


azimuth  =  it-  tan 


Pe 

Ps 
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elevation  =  tan  .  . 

Vp,’+p.’ 

Finally,  the  observ'ation  relation,  G ,  is  obtained  b>"  evaluating  equations  (3.18)  through  ( 3.23  ) 
sequentially,  followed  by  equation  (3.11). 

3.4.4  Observation  Matrix.  Once  the  reference  trajectory,  state  transition  matrix  and  observation 
relation  have  all  been  found,  the  next  step  is  to  determine  the  observation  matrix,  T .  The  observation 
matrix  does  for  the  observation  relation  what  the  state  transition  matrix  does  for  the  reference  trajectory. 

In  fact,  it  is  based  in  large  part  on  the  state  transition  matrix.  T  essentially  tells  the  algorithm  what 
change  in  the  predicted  observation  would  result  from  a  small  change  in  the  reference  state,  X .  The 
observation  matrix  is  defined  as 

T  =  H  O  ( 3.24 ) 


The  matrix  H  gives  the  change  in  the  observation  relation  accompanying  a  change  in  the  state  vector  at 
the  observation  time.  The  matrix  O  gives  the  change  in  the  state  vector  at  the  observ'ation  time 
accompanying  a  change  in  the  state  vector  at  the  reference  time.  Hence,  the  two  matrices  combine  to 
form  T ,  which  gives  the  change  in  the  observations  at  the  observation  time  caused  by  a  change  in  the 
reference  state  vector  at  the  reference  time. 


8G 

®*dX(t,) 

.  3X(t,) 

dx(0) 


T  =  HO  = 


dG  ax(ti) 

ax(tj)  ax(o) 


aG 

ax(o) 


(3.25) 


(3.26) 


But  in  this  case  the  observation  relation  and  state  transition  matrix  both  contain  only  the  position 
components  of  the  state  vector,  so  the  relation  becomes 


aG  aR(ti)  aG 
aR(ti)  ax(o)  ax(o) 


( 3.27  ) 
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As  described  equations  ( 3 .25 )  and  ( 3.27 ),  the  H  matrix  is  taken  as  the  partial  derivative  of 
the  observation  relation  with  respect  to,  (in  the  Taylor  series  method),  the  position  vector.  Thus,  the  H 
matrix  is  a  three  three  matrix  consisting  of  the  following  partial  derivatives. 


ac, 

ac, 

ac, 

aR, 

aRj 

3R3 

dG, 

ac^ 

3G, 

3R, 

aRj 

aRj 

3G, 

aG3 

aOj 

aR, 

3R, 

3R, 

Evaluation  of  these  partial  derivatives  involves  the  use  of  equations  ( 3.11 ),  ( 3.17 ),  (3.18 ),  ( 3.22 )  and 
( 3.23  )  as  well  as  judicious  use  of  the  chain  rule.  Again,  this  lengthy  development  will  not  be  presented 
here.  It  is  contained  in  its  entirety  as  Appendix  D. 

3.4.5  Instrument  Covariance  Matrix.  The  final  development  for  the  use  of  the  non-linear  least 
squares  algorithm  is  the  instrument  covariance  matrix,  Q .  The  instrument  covariance  matrix,  a  diagonal 
matrix,  is  obtained  from  the  reported  accuracies  of  the  radar  sites.  The  variances,  (the  square  of  the 
instrument’s  standard  deviations),  are  used  as  the  diagonal  elements. 

Q  = 

The  standard  deviations  used  in  this  thesis  are  shown  on  the  next  page. 

=  100  meters  =  1.567856x10"’  DUs 
g,  ■  =  0.025  deg  =  4.363323 x  10"^  rad 

^elevation  ^arinmlh 

3.4.6  Convergence  Criteria.  The  final  step  in  the  af^lication  of  the  least  squares  algorithm  is 
the  decision  of  a  convergence  criteria.  For  the  Taylor  series  method  of  track  compression  (and  throughout 
this  thesis)  the  correction  a^^lied  to  the  reference  state  is  compared  to  the  diagonal  elements  of  the  state 
covariance  matrix,  P ,  after  each  iteration.  If  the  correction  to  any  element  of  the  reference  state  is  larger 


range 

0 

0 


0 

_  2 
^  azimodi 

0 


0 
0 

^elevation  J 


(  3.29  ) 
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than  one  percent  of  the  square  root  of  the  corresponding  covariance  element,  another  iteration  is 
accomplished.  [20] 

3,4. 7  Summary.  To  summarize  the  detailed  methodology  of  the  Taylor  series  method  of  track 
compression,  a  track  of  data  is  split  into  arcs  of  no  more  than  90  seconds  in  length.  These  arcs  are  used  as 
observation  data  for  the  non-linear  least  squares  algorithm  described  in  chapter  two.  This  algorithm 
requires  the  development  of  a  reference  trajectoiy,  dynamics  model,  state  transition  matrix,  observation 
relation,  observation  matrix  and  instrument  covariance  matrix  as  well  as  a  decision  of  the  convergence 
criteria.  It  is  the  dynamics  model  that  primarily  distinguishes  the  Taylor  series  method  from  the 
integrator  method  to  be  described  in  chapter  four. 

The  non-linear  least  squares  algorithm  is  iterated  until  the  convergence  criteria  is  achieved 
Once  the  algorithm  has  converged  to  a  solution,  the  residuals  are  checked  to  ensure  that  they  are 
reasonable  (see  chapter  two).  The  algorithm  produces  a  reference  state  vector,  which  represents  the  state 
of  the  satellite  at  the  center  of  the  arc  of  data,  and  a  covariance  matrix,  which  represents  the  accuracy  of 
the  computed  reference  state. 


3.5  Results 

This  section  presents  the  results  of  the  Taylor  series  method  of  track  compression  portion  of  the 
thesis.  The  results  included  are  a  sample  of  the  program  output  and  the  first  and  last  pass  residuals  from 
track  compressions  of  several  different  orbits  and  viewing  geometries.  The  orbital  elements  used  for  the 
first  run  are  shown  in  Table  3-3. 


Table  3-3  Orbit  Nmnber  One 


semi-major  axis  (DUs) 

1.029994 

eccentricity 

0.003963 

inclination  (deg) 

40.000302 

1  longitude  of  ascending  node  (deg) 

19.998752 

argument  of  perigee  (deg) 

358.488511 

true  anomaly  (deg) 

46.513560 
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Radar  site  LION  was  used  for  this  first  run  with  observations  taken  every  0.4  seconds.  The 
satellite  was  visible  for  4  minutes,  37  seconds  producing  a  track  of  687  observations  which  were  split  into 
four  arcs.  Analyzing  the  residuals  shows  that  each  arc  was  successfully  compressed.  Figures  3-9  through 
3-11  show  the  residuals  from  the  first  pass  for  each  type  of  data,  for  the  first  of  the  foiu-  arcs.  These  are 
normalized  by  dividing  Ity  the  standard  deviations  for  each  type  of  data.  The  obvious  trends  in  the  data 
are  clear  indications  that  the  trajectory  evaluated  during  the  first  pass  is  not  the  optimal  trajectory.  The 
residuals  of  the  optimal,  converged  trajectory  should  appear  as  a  random  series  of  values  with  an  average 
near  zero  and  a  variance  near  one.  The  variance  is  near  one  because  the  residuals  are  normalized 
Approximately  68%  of  the  residuals  should  be  within  one  standard  deviation  of  the  mean. 


Figure  3-9  First  Pass  Range  Residuals 
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Figure  3-10  First  Pass  Azimuth  Residuals 


Figure  3-1 1  First  Pass  Elevation  Residuals 


The  last  pass  residuals  are  shown  in  Figures  3-12  through  3-14.  They  show  the  converged  solution,  with 
residuals  statistically  spread  as  described  previously. 
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Figure  3-12  Last  Pass  Range  Residuals 
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Figure  3-13  Last  Pass  Azimuth  Residuals 
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Figure  3-14  Last  Pass  Elevation  Residuals 


Similar  results  appeared  for  the  residuals  of  the  other  three  arcs.  To  summarize  these  results, 
once  the  convergence  criteria  has  been  achieved,  the  algorithm  calculates  the  average  normalized  residual 
for  each  type  of  observation  data  and  the  average  normalized  squared  residual.  These  are  recorded  and 
can  be  used  to  verify  the  validity  of  the  converged  solution.  The  following  table  shows  these  parameters 
from  this  first  run.  As  expected,  the  average  residuals  for  all  three  types  of  data  for  all  four  arcs  are  near 
zero  and  the  average  squared  residuals  are  near  one. 


Table  3-4  Orbit  One  Track  Compression  Results 


range 

azimuth 

elevation 

arc  1  residuals 

.00023364 

.00000251 

arc  1  squared  residuals 

.82982743 

.91883676 

arc  2  residuals 

gQj|j|j|gyy|jU 

.00000180 

-.00016240  1 

arc  2  squared  residuals 

1.02058824 

.91098848 

arc  3  residuals 

-.00028598 

.00000111 

.00004678 

arc  3  squared  residuals 

.98764160 

1.01218600 

arc  4  residuals 

-.00028314 

.00000415 

.00075569 

arc  4  squared  residuals 

.92798399 

.86142314 

1.03658611 

3-26 


The  results  in  Table  3-4  indicate  a  successful  track  compression.  The  maximiun  elevation 
recorded  in  the  data  was  less  than  four  degrees.  This  indicates  a  very  low  pass  in  the  sl^.  The  second 
orbit  was  selected  to  pass  quite  high  in  the  sky  over  the  COOK  radar  site.  Its  maximum  elevation  is  over 
68°.  The  following  shows  the  second  orbit  used  in  the  Taylor  series  method  of  track  compression.  It  also 
used  a  data  rate  of  one  observation  every  0.4  seconds. 


Table  3-5  Orbit  Number  Two 


semi-major  axis  (DUs) 

1.030011 

eccentricity 

0.003963 

inclination  (deg) 

39.997386 

longitude  of  ascending  node  (deg) 

225.000000 

argument  of  perigee  (deg) 

358.169601 

true  anomaly  (deg) 

46.835123 

This  orbit  is  quite  similar  to  orbit  one  but  with  the  plane  rotated  to  allow  it  to  pass  over  COOK. 
This  7  minute,  28  second  track  of  1 1 1 1  observations  was  split  into  five  arcs.  Table  3-6  shows  the  results 
of  the  track  compression  of  these  five  arcs. 

Table  3-6  Orbit  Two  Track  Compression  Results 


range 

azimuth 

elevation 

arc  1  residuals 

-.00003263 

-.00000854 

-.00012728 

arc  1  squared  residuals 

.88131956 

.95627818 

.97985336 

arc  2  residuals 

.00102424 

.00001146 

-.00019216 

arc  2  squared  residuals 

1.03146837 

.85304246 

.99571694 

arc  3  residuals 

-.00343549 

.00000031 

.01084358 

arc  3  squared  residuals 

.82358193 

.96928305 

.99036295 

arc  4  residuals 

.00100726 

-.00000616 

arc  4  squared  residuals 

.97109372 

1.00082980 

1.03648939 

arc  5  residuals 

-.00002229 

-.00000818 

-.00020686 

arc  5  squared  residuals 

.95489164 

.87828537 

.87487937 
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This  table  indicates  that  the  high  elevation  pass  also  converged  successfully.  The  next  step  was 
to  apply  the  Taylor  series  method  of  track  compression  to  a  higher  altitude  oibit.  For  this  case,  an  oibit 
two  Earth  radii  in  size  was  propagated  with  COOK  taking  observations  at  the  same  data  rate.  The 
satellite  passed  overhead  with  more  than  82°  of  elevation. 


Table  3-7  Oibit  Niunber  Three 


semi-major  axis  (DUs) 

2.029923 

eccentricity 

0.003921 

inclination  (deg) 

39.997902 

longitude  of  ascending  node  (deg) 

225.000684 

argument  of  perigee  (deg) 

357.730131 

true  anomaly  (deg) 

47.272337 

The  pass  lasted  94  minutes  producing  a  track  of  13,968  observations  split  into  63  arcs.  All  arcs 
were  successfully  compressed  with  valid  residuals.  All  of  these  residuals  were  near  the  desired  values  of 
one  and  zero.  T^le  3-8  shows  the  maximum,  average  and  minimum  of  the  63  average  residual  and 
average  squared  residual  results  for  the  range,  azimuth  and  elevation  of  the  63  arcs. 


Table  3-8  Orbit  Three  Track  Compression  Results 


maximiun 

average 

minimum 

average  range  residual 

0.00003117 

-0.00011385 

-0.00291641 

average  squared  range  residual 

1.17498481 

0.97460627 

0.76891445 

average  azimuth  residual 

0.00005437 

0.00000135 

-0.00004202 

average  squared  azimuth  residual 

1.20066433 

0.97300310 

0.78152302 

average  elevation  residual 

0.00123700 

-0.00018000 

-0.00898000 

average  squared  elevation  residual 

1.19814258 

0.98665555 

0.75231942 

None  of  these  values  are  too  far  from  the  nominal  values  of  zero  and  one.  The  final  orbit  was  eccentric 
with  the  same  semi-major  axis.  Again,  COOK  recorded  observations  at  the  same  data  rate. 
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Table  3-9  Orbit  Niunber  Four 


semi-major  axis  (DUs) 

2.030040 

eccentricity 

0.400019 

inclination  (deg) 

39.997573 

longitude  of  ascending  node  (deg) 

225.001504 

argument  of  perigee  (deg) 

IgiHglin 

1  true  anomaly  (deg) 

47.270383  | 

This  pass  lasted  2  hours  and  23  minutes.  The  track  consisted  of  21,296  observations  split  into  96 
arcs.  Table  3-10  presents  the  residual  results  in  the  same  manner  as  for  the  previous  orbit.  Again,  the 
values  are  all  reasonable  and  indicate  successful  convergence  on  all  arcs. 


Table  3-10  Orbit  Four  Track  Compression  Results 


1  II  maximum 

average 

minimum 

average  range  residual 

0.00004056 

-0.00029802 

-0.01989010 

average  squared  range  residual 

1-21929120 

0.98029666 

0.76489698 

average  azimuth  residual 

0.00011469 

0.00000372 

-0.00004421 

average  squared  azimuth  residual 

1.16064348 

0.96963272 

0.76134938 

average  elevation  residual 

0.00452605 

0.00004444 

-0.00474431 

average  squared  elevation  residual 

imggjQiyi 

3.6  Summary 

To  summarize  the  results,  the  Taylor  series  method  of  track  compression  is  quite  successful  at 
reducing  the  large  tracks  of  data  into  manageable  state  vectors  and  covariance  matrices.  The  drawback  to 
this  method  is  the  large  number  of  arcs  that  are  required  for  the  higher  altitude,  longer  track  passes.  The 
large  number  requires  more  time  to  compress  and  more  data  to  pass  to  the  SCC,  undermining  the  original 
intent  of  the  thesis.  The  next  portion  of  the  thesis,  integrator  track  compression,  was  designed  to 
eliminate  this  problem  1^  allowing  a  track  to  be  compressed  in  its  entirety,  without  being  split  into  arcs. 
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IV.  Integrator  Track  Compression 


4.1  Introduction 

This  chapter  details  how  the  integrator  method  of  track  compression  was  developed  using  the 
method  of  least  squares.  It  then  presents  the  results  of  the  integrator  track  compression  portion  of  the 
thesis.  The  notation,  units,  reference  fiames  and  truth  model  are  the  same  as  that  described  in  the  Taylor 
series  track  compression  description  in  the  previous  chapter. 

4.2  Detailed  Integrator  Track  Compression  Development 

The  integrator  method  of  track  compression  differs  from  the  Taylor  series  method  primarily  in 
the  dynamics  model.  This  model  numerically  integrates  the  equations  of  motion  and  the  equations  of 
variation  to  detennine  the  state  vector  and  state  transition  matrix  at  each  observation  time,  as  of^sed  to 
directly  calculating  these  as  in  the  T^Ior  series  method. 

4.2. 1  Comparison.  The  Taylor  series  method  was  originally  chosen  to  provide  for  faster 
computation.  The  suj^sed  advantage  dealt  with  the  necessity  of  a  numerical  integrator  to  evaluate  the 
state  frequently,  across  small  time  steps,  and  thus  prevent  a  direct  calculation  of  the  state  at  a  given 
observation  time.  However,  the  savings  in  time  did  not  manifest  itself  because  the  observation  times  are 
so  closely  spaced  that  a  numerical  integrator  can  go  directly  from  one  observation  time  to  the  next  without 
requiring  evaluation  across  the  interval.  This  fact,  coupled  with  the  rather  long  equations  that  developed 
out  of  the  Taylor  series  method,  allowed  the  numerical  integrator  to  produce  state  vectors  at  each 
observation  time  at  about  the  same  rate  as  the  Taylor  series  method  but  with  three  advantages.  The 
integrator  method  is  more  accurate,  is  easier  to  modify  to  account  for  additional  (fynamics  and,  most 
importantly,  does  not  require  the  track  to  be  split  into  arcs.  An  entire  track  can  be  compressed  at  once, 
and  only  one  state  vector  and  covariance  matrix  sent  to  the  SCC. 

4.2.2  Initial  Reference  State.  The  initial  reference  state  for  the  integrator  method  is  computed  in 
a  manner  very  similar  to  the  Taylor  series  reference  state.  The  reference  state  is  first  calculated  using  the 
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arc  and  line  method  described  in  chapter  three.  The  state  is  then  propagated  backward  to  the  time  of  the 
first  observation  in  the  track.  The  preliminary  orbit  determination  is  considered  more  accurate  in  the 
middle  of  the  track,  however,  the  integrator  method  lends  itself  to  the  use  of  a  reference  state  at  the 
beginning  of  the  track.  Therefore,  the  reference  state  is  calculated  at  the  middle  of  the  track,  using  the 
more  accurate  data,  and  then  propagated  to  the  beginning  of  the  track,  to  facilitate  the  integrator  (fynamics 
model. 

4.2.3  Dynamics  Model.  The  cfynamics  model  is  the  primary  difference  between  the  two  methods 
of  track  compression.  The  integrator  method  uses  the  same  fourth-order  Runge-Kutta  integrator  as  was 
described  for  the  truth  model  in  chapter  three. 

4.2.3. 1  Reference  Trajectory.  The  reference  trajectory  is  computed  by  propagating  the 
reference  state  through  each  observation  time.  This  propagation  requires  a  set  of  equations  of  motion 
describing  the  satellite’s  trajectory.  These  equations  take  into  account  two-bo^  acceleration  and  the 
effects  of  the  h  geopotential  harmonic.  They  are  the  same  as  was  used  for  the  truth  model. 

4.2.3.2  State  Transition  Matrix.  The  state  transition  matrix,  0(t,  ) ,  is  a  six  by  six 

matrix  used  to  describe  the  effects  of  a  small  change  in  the  reference  state  on  the  reference  trajectory.  As 
such,  it  is  different  for  each  observation  time.  Initially,  it  is  the  identity  matrix,  I .  It  changes  as  the 
interval  between  the  reference  time  and  the  observation  time  changes.  The  state  transition  matrix  is 
computed  via  the  following  method.  [7,  20] 


^(to,to)  =  I 

(4.1) 

^0(t,to)  =  A(t)0(t,to) 

(4.2) 

where 

^(t)“  Vxg|x„(,) 

andg  comes  from  the  equations  of  motion  of  the  state  vector. 
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^X  =  g(X,t) 


Thus,  the  components  of  A(t)  are  evaluated  as  the  partial  derivatives  of  the  equations  of  motion  with 
respect  to  the  state  vector.  These  are  known  as  the  equations  of  variation.  Once  again,  the  development 
of  these  equations  becomes  quite  lengthy  and  is  presented  as  Appendix  E. 

The  equations  of  motion  and  the  equations  of  variation  are  integrated  together.  The  equations  of 
motion  form  a  six  component  vector  of  equations  while  the  equations  of  variation  form  a  six  by  six  matrix 
of  equations.  Thus  the  propagator  must  integrate  42  equations  at  each  time  step.  This  seems  quite 
lengthy,  but  the  equations  are  not  as  long  as  the  Taylor  series  equations  and  this  propagation  provides 
both  the  reference  trajectory  and  the  state  transition  matrix.  Thus,  this  method  is  preferable  to  the  Taylor 
series  method. 

4.2.4  Taylor  Series  Similarities.  The  determination  of  the  observation  relation,  the  observation 
matrix  and  the  instrument  covariance  matrix  is  identical  to  the  description  contained  in  the  Taylor  series 
development.  The  convergence  criteria  is  also  identical. 

4.2.5  Summary.  To  summarize,  the  integrator  method  is  quite  similar  to  the  Taylor  series 
method.  The  difference  lies  in  the  (fynamics  model  where  a  fourth-order  Runge-Kutta  integrator  is  used 
to  numerically  integrate  the  equations  of  motion  and  the  equations  of  variation  to  provide  the  state  vector 
and  the  state  transition  matrix  at  each  observation  time.  This  results  in  a  more  accurate  state  vector 
derived  from  the  entire  track  of  observations. 

4.3  Results 

This  section  presents  the  results  of  the  integrator  method  of  track  compression  portion  of  the 
thesis.  The  results  will  be  presented  in  a  manner  similar  to  the  previous  chapter.  The  same  four  orbits 
were  used  to  test  the  integrator  track  compression  method  as  were  used  with  the  Taylor  series  method. 

The  first  set  of  orbital  elements  are  reprinted  as  Table  4-1. 
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Table  4-1  Orbit  Number  One 


semi-major  axis  (DUs) 

1.029994 

eccentricity 

0.003963 

inclination  (deg) 

40.000302 

longitude  of  ascending  node  (deg) 

19.998752 

argument  of  perigee  (deg) 

358.488511 

true  anomaly  (deg) 

46.513560 

The  entire  track  was  compressed  into  a  single  state  vector  and  covariance  matrix.  The  first  pass 
normalized  residuals  are  shown  in  Figures  4-1  tluough  4-3.  Just  as  in  the  last  chapter,  the  trends  indicate 
the  first  pass  was  not  the  optimal  solution.  All  three  sets  of  residuals  are  closer  to  zero  at  the  middle  of 
the  track.  This  occurs  because  the  reference  trajectory  used  to  generate  the  first  set  of  residuals  is 
constructed  fix>m  data  fi'om  the  middle  of  the  track. 


Figure  4-1  First  Pass  Range  Residuals 
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standard  deviations  standard  deviations 


0  100  200  300  400  500  600 

Figure  4-2  First  Pass  Azimuth  Residuals 


Figure  4-3  First  Pass  Elevation  Residuals 


Figures  4-4  throu^  4-6  show  the  last  pass,  converged  residuals.  As  expeaed,  they  average  near  zero  and 
aRiroximately  68%  are  within  one  standard  deviation  of  zero. 
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standard  deviations  standard  deviations 
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Figure  4-6  Last  Pass  Elevation  Residuals 


These  figures  indicate  the  compression  was  successful.  The  following  table  shows  the  average  residuals 
and  average  squared  residuals  (normalized)  for  each  data  Q^pe  for  the  first  orbit. 


Table  4-2  Orbit  One  Track  Compression  Results 


range 

azimuth 

elevation 

residuals 

.00056093 

,00001764 

.00367693 

squared  residuals 

.94139539 

.92858585 

1.00877648 

The  following  tables  complete  the  results  for  the  integrator  method  of  track  compression  portion 
of  the  thesis.  The  tables  present  the  last  three  orbits  (the  same  orbits  as  the  last  chapter)  and  the  residual 
data.  They  indicate  that  the  integrator  method  of  track  compression  is  also  a  valid  technique  for  data 
reduction  at  the  sites.  The  four  cases  combine  to  provide  a  range  of  orbit  types  and  viewing  angles. 
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Table  4-3  Orbit  Niunber  Two 


semi-major  axis  (DUs) 

1.030011 

eccentricity 

0.003963 

inclination  (deg) 

39.997386 

longitude  of  ascending  node  (deg) 

225.000000 

argument  of  perigee  (deg) 

358.169601 

true  anomaly  (deg) 

46.835123 

Table  4-4  Orbit  Two  Track  Compression  Results 


range 

azimuth 

elevation 

residuals 

-.00651527 

.00008277 

.02216073 

squared  residuals 

.94371461 

.93707563 

.98504838 

Table  4-5  Orbit  Number  Three 


semi-major  axis  (DUs) 

2.029923 

eccentricity 

0.003921 

inclination  (deg) 

39.997902 

longitude  of  ascending  node  (deg) 

225.000684 

I  argument  of  perigee  (deg) 

357.730131 

1  true  anomaly  (deg) 

47.272337 

. 

Table  4-6  Orbit  Three  Track  Compression  Results 


range 

azimuth 

elevation 

residuals 

-.00122617 

.00206993 

.03415182 

squared  residuals 

.98469824 

.98318182 

.99565793 
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Table  4-7  Orbit  Number  Four 


semi-major  axis  (DUs) 

2.030040 

eccentricity 

0.400019 

inclination  (deg) 

longitude  of  ascending  node  (deg) 

225.001504 

argument  of  perigee  (deg) 

357.731441 

true  anomaly  (deg) 

47.270383 

Table  4-8  Orbit  Four  Track  Compression  Results 


range 

azimuth 

elevation 

residuals 

-.04312409 

.02069728 

.12373816 

squared  residuals 

1.02945441 

1.00267420 

1.02216535 

4.4  Summary 

The  previous  results  indicate  that  the  integrator  method  of  track  compression  is  as  effective  as  the 
Taylor  series  method.  Additionally,  the  integrator  method  results  in  one  state  vector  and  covariance 
matrix  for  a  track.  This  reduces  the  data  recpiired  to  be  sent  to  and  processed  fay  the  SCC.  Another 
advantage  to  the  integrator  method  is  the  ease  with  which  the  djynamics  model  can  be  improved. 

Changing  the  equations  of  motion  is  a  relatively  simple  tasL  As  an  example,  additional  geopotential 
coefficients  and  atmospheric  drag  could  be  added. 

Although  the  results  of  the  track  compression  of  only  four  oibits  are  shown  here,  many  more 
were  actually  evaluated  for  this  thesis  (for  both  the  integrator  and  Taylor  series  methods).  The  final 
portion  of  the  thesis  involved  thousands  of  orbits  and  hundreds  of  thousands  of  track  compressions.  About 
one  track  per  ten  thousand  produced  compression  difficulties.  This  indicates  the  consistent  performance 
of  the  integrator  method. 
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y.  Global  Estimate 


5. 1  Introduction 

After  completion  of  the  two  methods  of  track  compression,  an  obvious  question  arose.  Can  the 
state  vectors  and  associated  covariance  matrices  from  multiple  compressed  tracks  be  used  to  provide  for  an 
accurate  orbit  estimate?  This  would  need  to  be  possible  if  the  track  compression  method  is  to  be 
implemented  by  the  SCC.  To  answer  this  question,  the  scope  of  the  thesis  was  extended.  This  new  effort 
was  called  Global  Estimate.  This  effort  involved  using  the  compressed  results  from  multiple  tracks  fi-om 
multiple  sites  around  the  world  as  input  data  to  another  least  squares  differential  corrector.  The 
application  of  the  least  squares  algorithm  differed  from  that  of  the  integrator  track  compression  method  in 
several  key  areas  to  be  discussed  in  the  following  sections. 

The  overall  goal  of  the  Global  Estimate  effort  was  to  demonstrate  that,  if  the  radar  sites  of  the 
SSN  forwarded  data  to  the  SCC  in  the  form  of  state  vectors  with  covariance  matrices,  instead  of  range, 
azimuth  and  elevation  readings  with  standard  deviations,  the  SCC  could  generate  accurate  orbit  estimates 
in  a  reasonable  period  of  time.  These  estimates  would  have  to  include  ballistic  coefficient  information. 

This  chapter  details  the  development  of  the  truth  model  used  to  generate  the  input  observation 
data  for  the  Global  Estimate  portion  of  the  thesis.  It  then  details  how  the  Global  Estimate  process  was 
developed  using  the  method  of  least  squares.  The  chapter  concludes  Ity  presenting  the  Global  Estimate 
results. 

Canonical  units  were  used  for  the  majority  of  the  Global  Estimate  development.  Exceptions  to 
the  canonical  units  standard  occur  in  the  time  steps  used  by  the  numerical  integrator  and  reporting  of 
track  compression  reference  times,  as  in  previous  sections,  as  well  as  the  ballistic  coefficient  and  air 
density.  The  ballistic  coefficient,  denoted  B*,  is  measured  in  units  of  m^/kg,  while  air  density  is 
calculated  in  imits  of  kg/m^. 
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5. 2  Detai led  Global  Estimate  Development 

This  section  describes  the  detailed  development  of  the  equations  and  matrices  used  for  the  Global 
Estimate  portion  of  the  thesis.  It  details  the  truth  model,  initial  reference  state,  and  estimator 
development. 

5. 2. 1  Truth  Model.  The  truth  model  used  for  the  Global  Estimate  portion  of  the  thesis  is  nearly 
identical  to  that  used  for  the  track  compression  portions.  The  difference  is  the  inclusion  of  air  drag  using 
an  atmospheric  model  presented  by  Regan.  [15]  The  truth  model  was  used  to  generate  the  radar  site 
observation  data  in  essentially  the  same  manner  as  for  the  two  track  compression  methods.  The  integrator 
method  of  track  compression  was  then  used  to  compress  the  observation  tracks  into  state  vectors  and 
associated  covariance  matrices.  These  did  not  include  ballistic  coefBcient  information.  The  equations  of 
motion  for  the  truth  model  Runge-Kutta  integrator  included  drag.  The  equations  of  motion  for  the  Global 
Estimate  dynamics  model  also  included  drag,  however,  the  integrator  in  the  track  compression  t^mamics 
model  did  not. 

5.2. 1. 1  State  Vector.  For  the  Global  Estimate  portion,  the  state  vector  is  extended  to 
allow  for  the  drag  calculations.  The  state  includes  a  term  representing  half  of  the  ballistic  coefficient. 

The  value  is  only  half  of  the  ballistic  coefficient  simply  to  save  the  computer  time  involved  in  calculating 
the  acceleration  due  to  drag,  which  will  be  shown  in  the  following  section  to  include  a  fector  of  Yi.  Thus 
the  seven-component  state  vector  for  these  portions  is  defined  as 

X  =  [r,  R,  R,;  V,  V,  V,.  (5.1) 

5. 2. 1. 2  Equations  of  Motion.  The  equations  of  motion  inclucte  the  acceleration  due  to 
the  geopotential  as  developed  in  chapter  three.  Additionally,  th^  include  the  acceleration  due  to  air  drag 
based  on  the  following  equation  [9] 

D  =  (5.2) 

where  D  is  the  force  due  to  drag 
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Cd  is  a  dimensionless  drag  coefficient  between  one  and  two  [16] 
p  is  the  local  atmospheric  density 

Vr  is  the  velocity  of  the  spacecraft  relative  to  the  local  atmosphere 
Af  is  the  exposed,  frontal  cross-sectional  area 

By  dividing  this  force  the  mass  of  the  satellite,  m,  the  acceleration  due  to  air  drag  is  shown  as 


R 


2  m 


(5.3) 


The  ballistic  coefficient  is  now  introduced  as  [2 1] 


B*  = 


m 


(5.4) 


The  velocity  relative  to  the  local  atmosphere  is  determined  by  assuming  an  atmospheric  model  which 
rotates  with  the  Earth.  Thus  the  relative  velocity  is  the  difference  between  the  satellite’s  inertial  velocity 
vector  and  the  cross  product  of  the  Earth’s  rotational  vector  with  the  satellite’s  position  vector.  By 
combining  this  with  equations  ( 5.3 )  and  ( 5.4 ),  and  expressing  the  result  in  vector  form,  the  following 
equation  for  the  acceleration  of  the  satellite  due  to  air  drag  is  obtained 

R  =  -^B*p  |v-(<d^xR)|  [V-{a)eXR)]  (5.5) 

where  <0®  =  0.05883359980154919  k  [22]. 

The  satellite-dependent  ballistic  coefficient  is  constant  and  is  included  as  part  of  the  state  vector. 
Specifically, 

X(7)sB^  (5.6) 

and 

Thus,  the  equations  of  motion  are  as  shown  on  the  next  page. 


« 
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Vi 

Vx 


x  = 


^gl  ■*"  ^dl 


^gK  •^dK 

0 


(5.7) 


where  Ag  is  the  acceleration  due  to  gravity  given  Ity  equation  ( 3.8 )  and  Ad  is  the  acceleration  due  to  air 
drag  given  by  equation  ( 5.5  ).  The  local  atmospheric  density  is  determined  using  an  atmospheric  model 
described  in  Aj^ndix  A.  [15] 

5. 2. 1.3  Observation  Data.  The  same  radar  sites  described  in  chapter  three  were  used 
for  this  portion  of  the  thesis.  For  this  case,  however,  all  sites  tracked  the  satellite  whenever  it  was  in  view. 
The  sites  produced  observations  at  the  rate  of  two  per  second.  The  tracks  were  then  compressed  using  the 
integrator  method  of  track  compression  described  in  chapter  four,  and  the  results  (state  vectors  and 
associated  covariance  matrices)  were  used  as  the  input  observation  data  for  the  Global  Estimate  least 
squares  estimator. 

5.2.2  Initial  Reference  State.  To  start  the  least  squares  algorithm,  an  initial  reference  state  is 
required.  For  Global  Estimate,  the  first  six  elements  of  the  initial  reference  state  are  taken  as  the  state 
vector  from  the  first  track  compression.  The  seventh  element  is  arbitrarily  assigned  a  value  of  0.05.  After 
one  iteration,  the  ballistic  coefficient  will  typically  be  corrected  to  a  value  with  at  least  two  digits  in 

with  the  final,  converged  value.  Thus,  the  aibitrary  nature  of  this  initial  value  does  not  hinder 
the  algorithm  and  convergence  is  still  rapidly  achieved. 

5.2.3  Dynamics  Model.  The  (tynamics  model  used  for  the  estimator  is  identical  to  that  used  by 
the  truth  model  integrator.  The  EOMs  are  the  same  as  equation  ( 5.7 )  described  in  section  5.2. 1.2.  The 
development  of  the  state  transition  matrix  requires  the  evaluation  of  the  A(t)  matrix.  This  requires  a 


development  similar  to  that  given  as  Appendix  E.  The  addition  of  the  drag  parameter  and  the 
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atmospheric  model  complicate  this  development.  The  resulting  equations  of  variation  are  developed  in 
their  entirety  in  Appendix  F. 

As  in  the  integrator  track  compression  method,  the  equations  of  motion  and  the  equations  of 
variation  are  integrated  simultaneously.  Because  of  the  additional  component  of  the  state  vector,  these 
novr  constitute  56  first-order  differential  equations. 

5.2.4  Observation  Relation.  The  observation  relation  is  considerably  simpler  in  Global  Estimate 
than  in  the  previous  developments.  The  input  data  consists  of  state  vectors  containing  three  position 
elements  and  three  velocity  elements  which  conespond  exactly  to  the  first  six  components  of  the  Global 
Estimate  state  vector.  The  input  data  gives  no  direct  information  about  the  seventh  component,  the  drag 
parameter.  Thus,  the  observation  relation,  G ,  is  a  direct  relation  of  the  first  six  state  vector  components. 


X,(t) 

x,(t) 

x,(t) 

x,(t) 

X,(t) 

x.(t). 


(5.8) 


5.2.5  Observation  Matrix.  The  observation  matrix,  T ,  is  again  given  by  the  following  relation. 


T  =  HO  = 


8G  8X(ti) 

8X(ti)  8X(0) 


8G 

8X(0) 


(3.26) 


Now,  however,  the  H  matrix  is  much  simpler. 


1  0  0  0  0  0  0 
0  1  0  0  0  0  0 
0  0  1  0  0  0  0 
0  0  0  1  0  0  0 
0  0  0  0  1  0  0 
0  0  0  0  0  1  0 


(5.9) 


5.2.6  Instrument  Covariance  Matrix.  The  instrument  covariance  matrix  changes  for  each 
observation,  as  o{^sed  to  the  constant  Q  of  the  track  compression  sections.  The  track  compression 
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algorithms  produce  both  a  state  vector,  used  as  the  input  data  to  Global  Estimate,  and  a  covariance 
matrix,  P .  This  covariance  matrix  is  different  for  each  compressed  track  and  serves  as  the  instrument 
covariance  matrix,  Q ,  for  Global  Estimate. 

5.2. 7  Convergence  Criteria.  The  convergence  criteria  differs  from  that  of  the  previous  chapters 

only  in  the  number  of  elements  evaluated.  The  algorithm  iterates  until  each  element  of  the  state  vector  is 
corrected  less  than  one  percent  of  the  square  root  of  the  corresponding  diagonal  element  of  the 

covariance  matrix.  This  time,  this  includes  the  seventh  diagonal  element  which  corresponds  to  the 
ballistic  coefficient  parameter. 

5.2.8  Summary.  To  summarize,  the  G/oW  Esft/nate  effort  was  aimed  at  demonstrating  the 
useffilness  of  the  track  compression  techniques  to  the  overall  orbit  estimation  problem.  The  project 
simulated  the  combination  of  data  forwarded  to  the  SCC  by  many  sites  all  tracking  the  same  satellite.  The 
same  non-linear  least  squares  differential  correction  algorithm  was  used  for  Global  Estimate  as  was  used 
for  the  track  compression  techniques,  but  the  implementation  was  considerably  different.  Most  notably^ 
the  (fynamics  model  requires  propagating  56  equations  of  variation  through  several  orbits  at  each 
iteration,  rather  than  propagating  42  equatiorts  through  a  relatively  short  arc,  as  in  the  integrator  method 
of  track  compression.  The  Global  Estimate  development  produces  a  state  vector  which  includes  a  drag 
parameter,  equal  to  one-half  the  ballistic  coefficient.  This  state  vector  is  derived  from  observations  spaced 
across  the  satellite’s  orbit  and  is  therefore  considered  orbital  estimation. 

5.3  Results 

This  section  presents  the  results  of  the  Global  Estimate  portion  of  the  thesis.  Results  from  foin 
orbits  are  included.  The  first  orbit  is  nearly  circular  and  nearly  equatorial.  The  second  is  slightly 
eccentric  and  retrograde.  The  third  orbit  is  larger  and  more  eccentric.  The  fourth  orbit  is  very  low  with  a 
high  ballistic  coefficient,  resulting  in  large  drag  effects.  The  Global  Estimate  algorithm  produces  a  state 
vector  and  covariance  matrix  at  a  specified  reference  time.  In  all  cases,  the  propagator  reported  the  actual 
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state  vector  at  this  reference  time  for  comparison.  These  results  include  the  comparison  between  the 
Global  Estimate  prediction  of  the  satellite’s  state  and  the  actual  state. 

The  following  table  shows  the  elements  of  the  first  orbit. 

Table  5-1  Orbit  Number  One 


semi-major  axis  (DUs) 

1.088173 

eccentricity 

0.000914 

inclination  (deg) 

5.540376 

longitude  of  ascending  node  (deg) 

37.686167 

argument  of  perigee  (deg) 

69.042996 

true  anomaly  (deg) 

254.358992 

ballistic  coefficient  (m^/kg) 

0.033730 

This  orbit  was  propagated  for  8  hours,  17  minutes,  36  seconds.  During  this  time,  19  tracks  of 
observations  were  recorded  and  compressed.  After  compression,  the  Global  Estimate  algorithm 
successfully  estimated  the  satellite’s  orbit.  The  following  table  shows  the  predicted  and  actual  state 
vectors,  as  well  as  the  difference.  The  state  vectors  are  presented  in  canonical  units  (with  the  excei^on  of 
the  ballistic  coefficient).  The  errors  are  presented  in  meters,  centimeters  per  second  and  square  meters  per 
kilogram. 


Table  5-2  Orbit  One  State  Vectors 


predicted 

actual 

error  | 

X-Position  (DUs) 

0.302786245 

0-302786054 

1.218  m  1 

Y-Position  (DUs) 

1.041225660 

1.041225678 

-0.115  m 

Z-Position  (DUs) 

0.065592718 

2.430  m 

X-Velocity  (DUs/TU) 

-0.921111898 

-0.921112005 

0.085  cm/s 

Y-Velocity  (DUs/TU) 

0.261496070 

0.261495883 

0.148  cm/s 

Z- Velocity  (DUs/TU) 

0.072230072 

0.072230044 

0.022  cm/s 

B*  (m"/kg) 

0.040041509 

0.033729594 

0.006  mVkg 

This  shows  the  accuracy  with  which  Global  Estimate  is  able  to  correctly  determine  the  satellite’s  orbit. 
Additionally,  the  covariance  matrix  is  computed.  The  square  roots  of  the  diagonal  elements  of  the 
covariance  matrix  roughly  correspond  to  the  three  axes  of  the  error  ellipsoids.  Thus,  the  first  three 
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diagonals  can  be  used  to  determine  the  expected  position  error,  and  the  next  three  can  be  used  to 
determine  the  expected  velocity  error.  The  following  table  reports  these  expected  errors  along  with  the 
magnitudes  of  the  actual  errors. 


Table  5-3  Orbit  One  Covariance  Statistics 


error  ellipsoid 

actual  error 

position  (m) 

3.45  ^ 

2.72 

velocity  (cm/s) 

3.92 

0.17 

The  second  orbit  was  retrograde.  Table  5-4  gives  its  orbital  elements.  This  was  propagated  for  6 
hours,  19  minutes,  28  seconds.  Nineteen  tracks  were  recorded.  Tables  5-5  and  5-6  show  the  results. 


Table  5-4  Orbit  Number  Two 


semi-major  axis  (DUs) 

1.160395 

eccentricity 

0.016347 

inclination  (deg) 

117.685067 

longitude  of  ascending  node  (deg) 

232.270593 

argument  of  perigee  (deg) 

93.762641 

true  anomaly  (deg) 

96.317677 

ballistic  coefficient  (m^/kg) 

0.072933 

Table  5-5  Orbit  Two  State  Vectors 


I — 

predicted 

actual 

error 

1  X-Position  (DUs) 

-0.798729115 

-0.798729336 

1.410  m 

Y-Position  (DUs) 

-0.469392066 

-0.469392202 

0.867  m 

Z-Position  (DUs) 

0.676601563 

0.676601794 

-1.473  m 

X-Velocity  (DUs/TU) 

0.124985690 

0.124985946 

-0.202  cm/s 

Y-Velocity  (DUs/TU) 

0.699969881 

0.699969644 

0.187  cm/s 

Z-Velocity  (DUs/TU) 

0.612644198 

0.612643960 

0.188  cm/s 

B*  (m^/kg) 

0.077966342 

0.072932682 

0.005  mVkg 
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Table  5-6  Orbit  Two  Covariance  Statistics 


error  ellipsoid 

actual  error 

position  (m) 

3.54 

2.22 

velocity  (cm/s) 

3.83 

0.92 

These  tables  indicate  that  Global  Estimate  had  no  problem  determining  the  orbit  of  this  retrograde 
satellite. 


The  third  orbit  attempted  was  much  larger.  Its  elements  are  given  in  Table  5-7. 

Table  5-7  Orbit  Number  Three 


1  semi-major  axis  (DUs) 

2.000000 

1  eccentricity 

0.400000 

1  inclination  (deg) 

28.702842 

longitude  of  ascending  node  (deg) 

75.373285 

argument  of  perigee  (deg) 

183.806110 

true  anomaly  (deg) 

304.719542 

ballistic  coefficient  (m^/kg) 

0.077227 

Orbit  three  was  propagated  for  4  hoiu^  and  48  minutes  or  one-fifth  of  a  day.  There  were  ten  tracks  of  data 
recorded  and  compressed.  Tables  5-8  and  5-9  show  the  results. 

Table  5-8  Orbit  Three  State  Vectors 


predicted 

actual 

error 

X-Position  (DUs) 

1.479112063 

1.479112334 

-1.730  m 

Y-Position  (DUs) 

-0.208099795 

-0.208099705 

-0.574  m 

Z-Position  (DUs) 

-0.811761258 

-0.8117610r2 

-1.569  m 

X-Velodty  (DUs/TU) 

0.395009244 

0.395009444 

“0.158  cm/s 

Y-Velocity  (DUs/TU) 

0.713048516 

0.713048427 

0.070  cm/s 

Z-Velocity  (DUs/TU) 

-0.108714162 

-0.108713749 

-0.326  cm/s 

B*  (m^/kg) 

0.061626080 

0.077227224 

-0.016  mVkg 
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Table  5-9  Orbit  Three  Covariance  Statistics 


error  ellipsoid 

actual  error  j 

position  (m) 

4.63 

2.41  1 

velocity  (cm/s) 

2.58 

0.37 

These  tables  indicate  that  Global  Estimate  was  also  able  to  successfully  estimate  the  orbit  of  this  higher 
altitude  satellite. 

The  fourth  and  final  satellite  had  a  lower  perigee  and  a  higher  ballistic  coefficient.  These 
combine  to  produce  a  high  level  of  drag  on  the  satellite.  Table  5-10  shows  the  orbital  elements  for  orbit 
four. 


Table  5-10  Orbit  Number  Four 


semi-major  axis  (DUs) 

1.053684 

eccentricity 

0.019270 

inclination  (deg) 

-  17.205573 

longitude  of  ascending  node  ((teg) 

118.168059 

argument  of  perigee  (deg) 

250.330578 

true  anomaly  (deg) 

307.348639 

ballistic  coefficient  (m^/kg) 

0.128213 

Orbit  four  was  propagated  for  half  of  a  d^.  There  were  3 1  tracks  of  data.  Tables  5-11  and  5-12  give  the 
results. 


Table  5-1 1  Orbit  Four  State  Vectors 


predicted 

actual 

error 

X-Position  (DUs) 

-0.007779346 

-0.007779829 

3 . 081  m 

Y-Position  (DUs) 

-1.041248918 

-1.041248650 

1.709  m 

Z-Position  (DUs) 

0.133952657 

0.133952453 

1.301  m 

X-Velocity  (DUs/TU) 

0.938547552 

0.938547812 

-0.206  cm/s 

Y-Velocity  (DUs/TU) 

-0.027152032 

-0.027152411 

0.300  cm/s 

Z- Velocity  (DUs/TU) 

-0.261756518 

-0.261756653 

0.107  cm/s 

B*  (m^/kg) 

0..  142901590 

0.128213036 

0 . 015  m^/kg 
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Table  5-12  Oibit  Four  Covariance  Statistics 


error  ellipsoid 

actual  enor 

position  (m) 

3.47 

3.76 

velocity  (cm/s) 

3.43 

0.38 

These  tables  show  that  Global  Estimate  was  also  successful  with  the  high  drag  case. 


5.4  Summary 

To  summarize,  the  Global  Estimate  effort  produced  an  algorithm  which  takes  compressed  tracks 
of  data  and  combines  them  using  the  non-linear  least  squares  differential  correction  method  to  estimate  a 
satellite’s  orbit.  The  algorithm  was  successful  at  determining  the  orbits  of  four  case  satellites  shown  in 
this  chapter.  Additionally,  many  thousands  of  additional  orbits  were  evaluated  in  the  simulation  effort, 
but  these  results  are  not  included  here. 
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VI.  Simulation 


6.1  Introduction 

After  completion  of  the  Global  Estimate  project,  another  question  arose.  Could  the  same  method 
of  orbit  estimation  employed  as  Global  Estimate  be  applied  to  an  entire  catalog  of  satellites?  That  is, 
could  the  sites  of  the  SSN  observe  the  satellites  frequently  enough,  and  the  SCC  process  the  data  quickly 
enough  to  keep  current  an  entire  catalog  of  seven  to  eight  thousand  objects?  The  Global  Estimate  project 
modeled  every  site  in  the  net^vork  tracking  a  single  object  every  time  it  was  visible.  In  the  real  world, 
such  a  situation  would  not  exist.  Sites  have  a  limited  tracking  ability  so  not  all  satellites  will  be  tracked 
when  the  geometry  allows.  The  visible  satellites  must  be  prioritized. 

To  answer  these  questions,  the  scope  of  the  thesis  was  again  broadened  to  include  large-catalog 
simulations.  To  keep  the  simulations  tractable,  160  satellites  (approximately  one-fiftieth  of  AFSPC’s 
catalog)  were  modeled.  The  sites  of  the  SSN  can  Really  track  50  objects  simultaneously.  The 
simulations  modeled  sites  with  the  ability  to  track  only  one  satellite  at  a  time.  This  way,  satellites  would 
be  tracked  about  as  fiequently  as  if  a  full-size  catalog  was  modeled  with  sites  tracking  50  satellites 
simultaneously.  This  provides  realistic  orbit  estimation  accuracy. 

The  simulation  portion  of  the  thesis  was  broken  down  into  three  phases.  The  objective  of  the  first 
phase  was  to  verily  the  concept  of  a  large-scale  simulation  and  to  verify  the  computer  code  used.  Several 
different  sets  of  orbits  were  simulated  to  observe  the  response  of  the  orbital  estinration  process  in  different 
altitude  regions.  The  second  phase  was  aimed  at  increasing  the  accuracy  of  the  covariance  matrices 
through  the  use  of  a  fading  memory  filter.  Phase  three  investigated  the  orbital  estimation  perforrttance 
under  conditiorrs  of  unpredictable  solar  and  atmospheric  activify.  This  chapter  describes  the  development 
of  and  the  results  of  the  three  phases. 
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6.2  Phase  One 


The  simulations  of  phase  one  amounted  to  the  propagation  of  160  satellites  through  several  days. 
Radar  sites  tracked  the  satellites  according  to  their  priority,  and  compressed  the  observation  tracks  into 
state  vectors  with  covariance  matrices.  These  compressed  tracks  were  used  to  estimate  the  orbits  multiple 
(typically  five)  times  daity.  The  estimates  were  kept  in  a  catalog,  propagated  forward  in  time  and  used  for 
the  prioritization  of  the  next  tracking  period.  Each  of  the  major  portions  of  the  simulation  is  detailed  in 
the  following  sections. 

6.2. 1  Orbits  Modeled.  Oibital  elements  and  drag  parameters  were  ran<tomly  generated  to 
provide  a  mix  of  different  types  of  realistic  orbits.  The  orbits  were  typically  low-Earth,  nearly  circular 
orbits.  Equations  (6.1)  were  used  to  determine  the  elements.  The  term  random  signifies  a  randomly 
generated  number  between  zero  and  one.  The  term  Gaussian  signifies  a  randomly  generated  number 
obQing  a  Gaussian  distribution  with  an  average  value  of  zero  and  a  variance  of  one.  The  minimum 
radius  was  set  at  200  kilometers  to  prevent  the  generation  of  orbits  that  dipped  too  far  into  the 
atmosphere. 


Rr  = 


Abs(Gaussian) 

10 


+  \finimum  Radius 


(6.1) 


R2  = 


Abs(Gaussiaii) 

iio 


+  Minimum  Radius 


A  = 


R]  +  R2 


ecc  =  Abs 


V 


2A 


f  TT 

inc  =  Mod|^Abs(Gaussian)x— ,jc 


Q  =  random  x  2k 
CO  =  random  X  2  it 
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V  =  random  x  2ji 


B*  =  Abs  {0.055+0.05  Gaussian) 

where  A  is  the  semi-major  axis 
ecc  is  the  eccentricity 
inc  is  the  inclination 
Cl  is  the  longiUuk  of  the  ascending  node 
(0  is  the  argument  of  perigee 
V  is  the  true  anomaly 
B*  is  the  ballistic  coefficient 

These  formula  provide  for  a  range  of  semi-major  axes  and  eccentricities  which  do  not  allow  the  satellites’ 
orbits  to  intersect  the  Earth’s  surface.  They  provide  a  range  of  inclinations,  mostly  less  than  45°  but 
possibly  as  high  as  completely  retrograde.  The  next  three  elements  are  all  angles  randomly  distributed 
between  zero  and  360°.  Finally,  the  radius  and  ballistic  coefficient  formulas  were  varied  to  allow  the 
study  of  different  ranges  of  perturbation  effects,  but  always  retained  the  same  basic  form  as  above.  The 
elements  of  a  set  of  160  satellites  generated  from  equations  (  6.1 )  and  used  during  a  typical  simulation  are 
shown  in  Figures  6-1  through  6-9  and  listed  in  a  table  as  Ai^ndix  G. 


Figure  6-1  Semi-Major  Axis 
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Earth  Radii  Earth  Radii 
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Figure  6-4  Eccentricity 
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Figure  6-5  Inclination 
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Figure  6-7  Argument  of  Perigee 


Figure  6-8  True  Anomaly 


6.2.2  Propagation.  After  the  elements  were  randomly  generated,  the  corresponding  seven- 
element  state  vectors  were  determined  using  the  method  described  in  section  2.5  of  Bate,  Mueller  and 
White.  [1]  The  satellites’  orbits  were  then  propagated  through  one  time  period.  Usually,  this  period  was 
one-fifth  of  a  d^.  In  one  case,  the  length  was  doubled.  All  160  satellites  in  the  catalog  were  propagated 
simultaneously.  The  radar  sites  tracked  the  satellite  in  view  with  the  highest  priority.  The  prioritization 
scheme  will  be  described  in  section  6.2.8.  Observation  data  was  generated  using  the  method  described  in 
chapter  three. 

Additionally,  four  more  zonal  harmonics,  J3  -  Je,  were  included  in  the  truth  model.  The 
development  of  the  equations  of  motion  due  to  these  effects  proceeds  identically  to  that  of  J2  shown  in 
equations  ( 3.6 )  through  ( 3.8 ).  This  lengthy  development  is  provicted  in  Appendix  B. 
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6.2.3  Track  Compression.  Whenever  the  satellite  being  actively  tracked  Ity  a  site  set  (passed 
below  the  local  horizon  of  the  site)  the  track  compression  program  was  called  to  compress  the 
observations  into  the  same  type  of  state  vector  and  covariance  matrix  described  previously.  The  integrator 
method  of  track  compression  was  again  used  for  the  simulation. 

6.2.4  Global  Estimate.  After  the  propagator  had  generated  observation  data  for  one  period,  and 
the  track  compression  subroutine  had  compressed  all  of  this  data  into  state  vectors  and  covariance 
matrices,  the  Global  Estimate  method  of  orbit  estimation  was  applied  to  each  of  the  160  satellites.  The 
results  of  Global  Estimate  were  retained  as  the  catalog.  This  included  the  reference  time,  the  seven- 
element  state  vector  and  the  49 49  element  covariance  matrix. 

6.2.5  Bayes  Filter.  For  the  second  and  subsequent  estimations,  a  slightly  different  method  of 
orbit  ^termination  was  required.  The  method  chosen  is  the  Bayes  filter  described  in  chapter  two.  This  is 
identical  to  the  least  squares  method,  except  that  the  results  of  a  previous  estimate  are  used  as  both  the 
initial  reference  trajectory  and  as  an  observation  to  be  processed.  This  observation  differs  from  the  others 
only  in  the  level  of  accuracy  and  the  inclusion  of  the  drag  parameter.  The  Bayes  filter  produces  the  same 
type  of  state  vector  and  covariance  matrix  as  Global  Estimate. 

The  observation  matrix  for  the  Bayes  filter  is  different  when  using  the  previous  estimate  as  a 
piece  of  data.  In  this  case,  the  matrix  is  a  seven  by  seven  identity  matrix.  The  instrument  covariance 
matrix  now  becomes  the  seven  by  seven  covariance  matrix  of  the  previous  estimate.  The  state  vectors  and 
covariance  matrices  for  the  remainder  of  the  data  (from  the  newly  compressed  tracks)  are  handled  exactly 
as  described  in  the  Global  Estimate  method. 

6.2.6  Update  Satellite.  After  the  Global  Estimate  or  Bayes  filter  method  of  orbit  estimation  has 
been  aiq)lied  to  all  160  satellites,  the  catalog  is  updated  to  the  end  of  the  propagation  cycle.  This  is 
required  because  the  two  estimation  processes  provide  results  for  each  satellite  at  the  same  reference  time 
as  each  satellite’s  last  track  compression.  At  the  end  of  the  estimation  portion  of  the  cycle,  the  propagator 
is  at  the  end  of  the  cycle  time,  but  the  estimates  in  the  catalog  are  at  the  various  reference  times  of  the 
track  compressions.  The  update  satellite  portion  propagates  the  state  vector  and  covariance  matrix  for 
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each  satellite  from  its  catalog  reference  time  to  the  same  reference  time  as  the  simulator,  which  is  the  end 
of  the  one-fifth  day  period 

6.2. 7  Truth  Catalog  Comparison.  Once  the  catalog  of  estimates  has  been  propagated  to  the  end 
of  the  cycle,  it  is  compared  to  the  truth  data.  For  each  satellite,  each  element  of  the  vector  is  compared  to 
determine  the  true  error  in  the  prediction.  The  results  of  this  comparison  are  the  key  output  desired  from 
the  simulation. 

6.2.8  Prioritize.  At  the  completion  of  each  cycle,  the  satellites  are  reprioritized  Each 
covariance  matrix  is  propagated  to  the  end  of  the  next  cycle  to  determine  how  big  the  covariance  would 
become  if  no  estimates  occurred  during  the  next  cycle.  The  satellites  are  rank  ordered  based  primarily  on 
covariance  matrix  size,  but  also  predicted  satellite  reentry  is  taken  into  account  and  given  high  priority.  If 
a  satellite’s  orbit  has  not  been  successfully  estimated,  it  is  not  propagated  Instead,  that  satellite  is  also 
assigned  a  high  priority. 

6.2.9  Output.  The  simulator  produces  several  output  files  at  various  points  of  the  cycle.  A  log 
file  is  kept  which  keeps  track  of  when  each  process  is  started  and  stopped,  the  significant  results  of  some 
processes  and  error  messages  generated  during  execution.  At  the  end  of  the  propagation  phase,  the 
catalog  of  truth  data  is  recorded  At  the  end  of  the  estimation  phase,  the  catalog  of  estimates  is  recorded. 
At  the  end  of  each  cycle,  a  file  is  generated  recording  the  differences  between  the  true  state  vectors  and  the 
estimated  state  vectors.  These  differences  are  also  normalized  by  the  associated  covariance  element  and 
recorded.  These  files  combine  to  provide  a  detailed  record  of  the  simulation. 

6.2.10  Simulation  Robustness.  Several  pieces  of  code  were  specifically  modified  to  prevent 
program  crashes  during  long-term  simulations.  The  majority  of  these  modifications  were  designed  to 
prevent  the  attempted  propagation  of  a  satellite  through  regions  of  the  atmosphere  too  low  to  be  effectively 
described  by  the  atmospheric  model.  These  modifications  included  preventing  a  compressed  track  of  data 
from  being  used  as  the  initial  reference  trajectoiy  during  a  Global  Estimate  run  if  the  track  of  data  was  too 
short  and  hence  not  of  sufficient  accuracy.  Additionally,  ballistic  coefficients  were  checked  after  each 
iteration  to  ensure  they  were  not  so  far  out  of  range  that  the  satellite  would  reenter  or  escape  during  the 
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next  iteration.  After  convergence,  the  resulting  ballistic  coefBcients  were  again  checked  against  a  more 
stringent  window.  This  prevented  such  obvious  anomalies  as  negative  values.  Additionally,  cotte  was 
added  which  watched  for  an  iteration  that  was  diverging.  If  three  successive  corrections  to  the  state  vector 
each  increased  in  magnitude,  the  estimation  was  postponed  until  more  data  was  available  (the  next  cycle). 
Round-off  error  occasionally  causes  a  diagonal  element  of  a  covariance  matrix  to  become  negative. 
Theoretically,  this  could  never  haj^n,  and  the  elements  are  assumed  positive  at  several  points  in  the 
simulation  where  the  square  root  of  the  element  is  needed.  To  prevent  crashes,  the  diagonal  elements 
were  monitored  to  ensure  that  no  attempt  was  made  to  evaluate  the  square  root  of  a  negative  covariance 
element. 

By  adding  code  to  prevent  crashes,  and  by  recording  frequent  data  files  for  restarting  in  case  it 
did  crash,  the  simulation  code  was  made  quite  robust.  As  a  result,  simulations  can  be  run  spanning 
several  weeks,  with  no  loss  of  data. 

6. 3  Phase  One  Results 

The  objective  of  phase  one  was  to  verify  the  integration  of  the  various  algorithms  used  previously 
in  the  thesis.  Additionally,  research  objectives  included  evaluating  the  effect  of  varying  levels  of 
atmospheric  drag  on  the  ability  to  estimate  the  orbits  of  the  catalog,  and  the  use  of  different  periods  for  the 
propagate-estimate-update  cycle.  The  elements  for  the  first  simulation  were  generated  using  equations 
(6.1).  The  elements  generated  are  shown  in  Figures  6-1  through  6-9  and  listed  as  Appendix  G.  The 
simulation  was  run  for  two  days. 

The  kQ^  parameters  monitored  during  each  simulation  will  now  be  described.  First,  the  number 
of  satellites  whose  orbits  have  been  estimated  is  tracked.  This  parameter  is  reported  as  a  fraction  of  the 
size  of  the  current  catalog.  The  size  of  the  catalog  changes  with  time  as  some  satellites’  orbits  will  decay 
to  reentry.  These  satellites  are  then  removed  from  the  catalog  and  their  statistics  are  no  longer  relevant. 
The  number  of  satellites  removed  due  to  reentry  is  reported  as  a  percentage  of  the  original  size  of  the 
catalog.  The  number  of  satellites  whose  positions  are  predicted  within  one  kilometer  of  their  true 
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positions  is  tracked  This  parameter  is  reported  both  as  a  fraction  of  the  catalog,  and  as  a  fraction  of  the 
number  of  estimated  satellites.  Obviously,  after  every  satellite  in  the  catalog  has  been  estimated,  these  two 
descriptions  become  identical.  Finally,  a  covariance  accuracy  parameter  is  tracked.  The  covariance 
matrices  computed  at  each  cycle  should  statistically  match  the  actual  estimation  accuracies  according  to 
the  three-dimensional  Gaussian  probability  distributions.  This  states  that  approximately  19.9%  of  the 
actual  error  vectors  should  lie  within  the  one-sigma  error  ellipsoid  reported  as  the  diagonal  elements  of 
the  covariance  matrices.  The  number  of  position  error  vectors  that  lie  within  their  respective  one-sigma 
error  ellipsoids  is  reported  as  a  fraction  of  the  munber  of  estimated  orbits.  If  the  (fynamics  were  perfectly 
described  1^  all  of  the  estimation  algorithms  used  in  the  oibit  estimation  process,  then  we  would  expect 
that  this  parameter  would  average  around  19.9%.  However,  since  the  (fynamics  used  Ity  the  truth  model 
are  not  perfectly  modeled  1^  the  estimator  algorithms  (as  in  real  life)  the  actual  errors  may  be  much 
larger.  These  five  percentages  are  all  plotted  at  the  completion  of  each  cycle.  These  same  parameters  are 
used  to  show  the  results  of  most  of  the  simulations  in  the  thesis.  Figure  6-10  shows  the  results  of  the  first 
simulation. 


Figure  6-10  Simulation  One  Results 
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Over  the  two-day  period,  five  satellites  reentered,  all  satellites’  oibits  were  estimated,  and  at  one  point, 
nearly  80%  of  the  catalog  was  estimated  within  one  kilometer.  However,  the  accuracy  of  the  estimations 
decayed.  By  the  end  of  the  second  day  of  tracking,  this  parameter  was  down  to  60%  and  was  still  falling. 
Additionally,  the  covariance  matrices  computed  were  far  too  optimistic.  This  manifests  itself  in  the  decay 
of  the  Vo  of  estimated  within  one  sigma  parameter  to  zero.  Both  of  these  problems  are  caused 
differences  in  the  dynamics  of  the  truth  model  and  the  estimator  models. 

The  next  set  of  oibits  was  generated  to  observe  increased  drag  effects.  The  following  two  figures 
show  the  altitudes  of  perigee  and  ballistic  coefficients  used  in  simulation  two. 


Figure  6-11  Altitude  of  Perigee  —  Simulation  Two 


Figure  6-12  Ballistic  Coefficient  -  Simulation  Two 


6-11 


The  lower  altitudes  and  increased  ballistic  coefficients  were  expected  to  increase  the  number  of  satellites 
which  reentered  and  decrease  the  accuracies  of  the  reported  covariance  matrices.  Figure  6-13  shows  the 
actual  results  observed. 


Figure  6-13  Simulation  Two  Results 

As  expected,  more  satellites  reentered  during  the  two  day  period,  and  the  accuracy  of  the  covariance 
matrices  dropped  dramatically.  The  accuracies  of  the  estimates  did  not  change  substantially  over 
simulation  one,  despite  the  increased  drag.  This  is  likely  because  many  of  the  satellites  with  more  drag 
have  reentered  1^  the  end  of  the  second  day,  thus  they  are  not  included  in  the  %  within  1  km  parameters. 
Also,  the  lower  orbits  allow  more  frequent  observation  opportunities,  which  improves  the  estimations. 

The  next  set  of  orbits  was  generated  to  see  the  effects  of  the  estimation  process  for  higher  altitude 
satellites.  The  same  ballistic  coefficients  as  simulation  one  were  used.  These  same  ballistic  coefficients 
were  used  for  all  other  simulations  of  all  three  phases  as  well.  The  altitudes  of  perigee  are  shown  in 
Figure  6-14. 
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Figure  6-14  Altitude  of  Perigee  —  Simulation  Three 
This  simulation  was  expected  to  be  less  accinate  than  the  previous  two,  as  less  frequent  tracks  of 
observations  are  recorded.  A  decline  in  reentries  was  also  ejqiected. 


Figure  6-15  Simulation  Three  Results 

During  the  two  day  simulation,  no  satellites  reentered,  and  it  took  the  full  two  days  to  estimate  every  orbit. 
This  is  because  of  the  lower  frequency  of  observation  tracks.  Also,  the  average  accuracy  was  lower,  but 
not  significantly. 
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The  next  simulation  was  designed  to  show  the  effect  of  increased  cycle  period.  Instead  of  the 
propagate-estimate-update  cycle  running  over  one-fifth  day  periods,  it  was  extended  to  two-fifth  day 
periods.  All  other  parameters  were  the  same.  Figure  6-16  shows  the  results. 


Figure  6-16  Simulation  Four  Results 

The  parameter  trends  are  basically  the  same,  indicating  that  the  period  is  not  much  of  a  fector  in 
determining  the  accuracy  of  the  estimates  or  of  the  covariance  matrices. 

The  fifth  simulation,  the  final  simulation  of  phase  one,  modeled  orbits  substantially  higher  than 
simulations  three  and  four.  The  size  of  the  orbits  stretched  to  nearly  three  Earth  radii,  as  shown  by  the 
perigee  altitudes  in  Figure  6-17. 


Figure  6-17  Altitude  of  Perigee  -  Simulation  Five 
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Because  it  took  much  longer  to  estimate  a  substantial  portion  of  the  catalog,  the  simulation  was  run  for 
three  days  instead  of  the  usual  two.  Figure  6-18  shows  the  results. 


day 

Figure  6-18  Simulation  Five  Results 

After  three  days,  the  accuracy  of  the  catalog  seemed  to  have  peaked  and  began  a  decline,  so  the  simulation 
was  terminated  despite  having  only  estimated  101  satellites. 

The  first  phase  of  the  simulation  showed  that  the  different  software  components  worked  together 
as  expected.  The  various  parameters  behaved  as  expected.  Thus,  having  achieved  the  objectives  of  phase 
one,  the  next  phase  was  entered. 

6.4  Phase  Two 

The  objective  of  phase  two  was  to  increase  the  accuracy  of  the  covariance  matrices.  Because  the 
(fynamics  models  of  the  track  compression  and  orbit  estimation  algorithms  differ  from  that  of  the  truth 
model,  the  actual  orbit  estimations  are  not  as  accurate  as  reported.  This  was  seen  in  the  majority  of  the 
phase  one  simulations.  The  covariance  matrices  must  be  accurate  if  they  are  to  be  used  operationally  to 
provide  a  confidence  level  for  the  estimates.  As  described  in  chapter  two,  a  fading  memory  filter  is  one 
way  to  attack  this  problem.  A  fading  memory  filter  allows  the  estimator  to  place  less  weight  on  older  data 
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and  also  allows  the  update  portion  of  the  cycle  to  increase  the  size  of  the  covariance  matrices  with  time, 
thus  making  them  more  realistic.  This  requires  the  determination  of  a  set  of  coefficients  to  use  as  the 
diagonals  of  the  p-matrix.  Phase  two  of  the  simulation  was  devoted  to  determining  the  correct  values  to 
use  for  the  seven  coefficients. 

Assumptions  were  made  as  to  the  structure  of  the  p-matrix.  The  same  value  was  used  for  the 
first  six  coefficients  and  a  value  of  one  was  used  for  the  seventh.  This  amounts  to  full  remembrance  of  the 
ballistic  coefficient  information,  and  equal  remembrance  of  each  of  the  other  state  elements. 


6.5  Phase  Two  Results 

The  first  value  used  (for  the  first  six  elements)  was  0.93.  This  equates  to  a  data  half-life  of  about 
nine  and  one-half  days.  The  same  orbital  elements  and  ballistic  coefficients  used  for  simulation  one  were 
used  for  all  of  phase  two.  The  results  of  simulation  six  (with  p  =  0.93)  are  shown  in  Figure  6-19. 


Figure  6-19  Simulation  Six  Results 

After  1.4  days,  the  simulation  was  stopped.  The  figure  clearly  shows  that  the  value  of  0.93  is  too  large, 

I 

i.e.  a  lower  value  should  be  used  to  allow  the  estimator  to  reduce  the  weight  of  the  older  data  even  more. 


6-16 


— □ —  %  reentered 
— O— %  of  catabg  within  1  km 
—A—  %  of  estimated  w  ithin  1  km 
— X— %  estimated 
— X~  %  of  estimated  w  ithin  1  sigma 


^ - □ - □ - o - □ - C3 - 1 - 1 - 1 - ! 

1  2 

day 

6-21  Simulation  Eight  Results 
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The  value  of  0.80  was  still  too  high.  The  next  figure  shows  the  results  when  the  value  of  0.50  was  tried. 
This  extreme  case  corresponds  to  a  half-life  of  only  one  day. 


Figure  6-22  Simulation  Nine  Results 

After  simulation  nine  was  stopped,  a  new  approach  was  evidently  required.  The  covariance 
matrices  were  fer  too  optimistic,  even  with  an  extreme  fading  memoiy  filter.  Additionally,  the  estimates 
themselves  were  not  as  accurate  as  expected.  The  suspected  cause  of  both  of  these  problems  was  the 
additional  dynamics  of  the  truth  model.  The  estimator  should  be  able  to  maintain  a  high  degree  of 
accuracy  even  without  all  of  the  ^mamics  included  in  the  truth  model,  as  long  as  the  additional  d^mamics 
do  not  introduce  a  systematic  error.  The  truth  model  can  even  include  a  random  acceleration  to  model  a 
stochastic  system,  and  the  estimator  should  still  produce  good  estimates.  The  suspected  problem  in  this 
case  was  the  inclusion  of  the  additional  zonal  harmonics  in  the  truth  model.  Additional  harmonics  will 
cause  the  oibital  period  to  decrease  very  slightly.  For  this  estimator,  the  size  of  the  orbit  is  being 
determined  so  precisely  that  this  small  change  in  period  is  noticeable.  If  this  were  the  cause  of  the 
problems,  then  the  actual  errors  would  likely  lie  in  the  direction  of  the  velocity  vectors.  To  investigate 
this  possibility,  the  following  plots  were  generated  to  show  the  component  of  each  actual  error  vector, 
after  one  day  of  propagation,  resolved  along  three  different  axes.  The  first  figure  shows  the  fiaction  of 
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each  error  vector  which  lies  in  the  direction  of  the  corresponding  velocity  vector.  The  second  shows  the 
traction  resolved  along  the  radial  direction.  The  third  shows  the  out  of  plane  fi^ion.  These  were 
determined  by  taking  the  dot  product  of  each  error  vector  with  a  unit  vector  in  each  of  the  three  directions. 
Because  the  orbits  are  not  perfectly  circular,  the  three  directions  are  not  orthogonal.  Nevertheless,  the 
three  figures  show  a  very  clear  trend. 
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majority  were  caused  Ity  the  satellite’s  position  being  estimated  behind  its  actual  position.  This  indicates 
that  the  inclusion  of  the  additional  zonal  harmonics  introduced  a  systematic  error  into  the  estimation 
process.  To  prove  that  the  zonal  harmonics  were  indeed  the  source,  an  additional  simulation  was 
completed  with  the  ^-coefficient  set  to  one  (lull  remembrance  of  all  elements)  and  the  values  for  the 
additional  zonal  harmonics  (Ja-Je)  set  to  zero  in  the  truth  model.  Figure  6-26  shows  these  results. 


Figure  6-26  Simulation  Ten  Results 
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Based  on  Figure  6-26,  it  was  decided  tliat  the  estimator  was  accurate  enough  to  see  the  effect  of 
the  additional  harmonics.  These  harmonics  introduced  a  systematic  error.  This  raised  a  very  important 
point  in  the  effectiveness  of  this  entire  process  of  orbital  estimation.  All  sources  of  systematic  errors  must 
be  determined  and  included  in  the  estimator  dynamics  model.  For  phase  three,  the  zonal  harmonic 
coefficients  were  left  at  zero  in  the  truth  model  to  remove  their  efiect.  In  an  operational  enviromnent,  the 
estimator  model  should  include  a  much  more  accurate  geopotential  model.  In  fact,  any  source  of 
systematic,  as  oi^sed  to  random,  error  should  be  included.  These  are  ^ically  deterministic  in  nature 
and  their  inclusion  is  not  particularly  difficult.  For  example,  luni-solar  effects  might  become  significant 
in  an  operational  orbit  estimator.  These  effects  are  easily  predicted  and  would  be  included  in  the 
estimator’s  ^mamics  model. 

The  objective  of  phase  two  was  to  increase  the  accuracy  of  the  covariance  matrices.  This 
objective  was  achieved,  although  not  through  the  use  of  the  fading  memory  filter  as  expected.  Instead,  the 
(fynamics  of  the  truth  mo(tel  had  to  be  modified  to  more  closely  match  those  of  the  estimator.  This  is  not 
the  most  desirable  situation,  however,  the  effects  of  the  additional  zonal  harmonics  could  be  included  in 
both  models.  Since  these  effects  are  completely  deterministic,  they  would  not  significantly  change  the 
results. 

6.6  Phase  Three 

The  first  two  phases  used  an  atmospheric  model  which  was  identical  in  both  the  truth  model  and 
the  estimator  model.  In  reality,  nearly  all  orbital  perturbations  can  be  predicted  except  for  air  drag.  This 
is  because  the  atmo^heric  density  varies  quite  significantly  both  spatially  and  temporally.  Solar  heating, 
solar  flares,  geomagnetic  storms  and  local  atmospheric  weather  all  affect  the  local  density  at  any  given 
time.  These  effects  are  largely  unpredictable,  leading  to  inaccurate  atmospheric  density  models.  To 
model  real  world  effects,  phase  three  of  the  simulation  included  a  different  atmospheric  model  in  the 
Ramies  of  the  truth  model.  The  estimator  (tynamics  remained  the  same.  This  allows  for  a  more 
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realistic  simulation  of  the  actual  results  which  could  be  expected  from  the  use  of  the  orbital  estimation 
process  with  real  data. 

The  original  atmospheric  model  described  in  Appendix  A  produces  a  local  density  based  on  the 
altitude  and  the  sea-level  atmospheric  pressure.  By  varying  the  sea-level  pressure,  slightly  different 
densities  are  calculated  The  variation  of  this  pressure  is  the  basis  for  the  alterations  of  the  atmospheric 
model  of  phase  three.  Three  different  effects  were  modeled. 

An  investigation  of  more  accurate  models  led  to  altering  the  sea-level  pressme  ±  25%  to  allow 
day-night  density  variations  more  in  line  with  actual  density  profiles.  [10, 14]  This  first  effect  modeled 
the  atmospheric  bulge  caused  by  solar  heating  on  the  sunward  side  of  the  Earth.  This  was  included  in 
simulation  eleven  by  increasing  or  decreasing  the  sea-level  pressure  by  up  to  25%  depending  on  where  the 
satellite  was  relative  to  the  sun.  For  purposes  of  the  simulation,  it  only  matters  that  the  bulge  is  in  a 
consistent  area  relative  to  the  Geocentric-Inertial  reference  frame.  Therefore,  the  bulge  was  modeled  in 
the  positive  direction  of  the  I-axis.  The  amount  of  increase  or  decrease  was  based  on  the  angle  between 
the  satellite’s  instantaneous  position  vector  and  the  I-axis.  By  changing  the  cfynamics  of  the  truth  model 
to  provide  this  25%  variation,  the  simulation  was  made  more  realistic. 

Additionally,  because  the  density  fluctuates  because  of  other,  unpredictable  factors,  an  additional 
random  variation  was  added.  The  sea-level  pressure  was  changed  by  a  random  percentage  following  a 
Craussian  distribution.  The  standard  deviation  of  this  distribution  is  5%.  The  third  variation  was  not 
included  in  simulation  eleven  and  will  be  discussed  later. 

Finally,  because  the  more  realistic  atmospheric  model  was  employed,  an  attempt  was  also  made 
to  generate  a  catalog  of  satellites  with  altitudes  which  more  closely  resembled  those  of  the  actual  catalog. 
84%  of  the  objects  in  the  actual  catalog  are  over  800  kilometers  in  altitude.  [4]  The  equations  generating 
the  catalog  were  modified  to  provide  a  distribution  with  84%  of  the  semi-major  axes  over  800  kilometers. 
The  following  histogram  shows  the  new  distribution  of  semi-major  axes. 
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attitude  of  semi-m^or  axis  (km) 

Figure  6-27  Phase  Three  Altitude  Distribution 
This  distribution  was  used  for  all  simulations  of  phase  three. 

Finally,  a  nmnber  of  changes  to  the  estimation  algorithm  were  implemented  to  make  the 
estimations  more  accurate.  Among  these  changes  was  a  check  of  the  residuals  after  the  oibit  was 
estimated.  If  the  residuals  were  unrealistic,  the  new  solution  was  rejected.  Additionally,  the  number  of 
successive  estimation  attempt  failures,  for  any  reason,  for  each  satellite  was  tracked.  The  priority  of  a 
given  satellite  was  raised  if  the  algorithm  failed  to  estimate  the  orbit.  This  provides  more  tracks  of  data 
for  the  problem  satellite  during  the  next  estimation  attempt.  Finally,  if  the  Bayes  filter  failed  to  arrive  at 
an  accurate  estimate  with  more  than  30  new  tracks  of  data,  the  previous  estimate  was  dropped  and  the 
least  squares  batch  estimator  was  used  to  provide  a  new  estimate.  If  this  was  not  successful,  the  oldest 
tracks  of  data  would  be  dropped,  one  at  a  time,  until  the  least  squares  estimator  was  successful,  or  until 
only  30  tracks  of  data  remained.  This  prevents  outdated  data  fi'om  causing  failed  estimate  attempts  when 
sufficient  new  data  is  alreatfy  available  to  provide  an  accurate  oibit  determination. 

6. 7  Phase  Three  Results 

Figure  6-28  shows  the  results  of  simulation  eleven.  Again,  simulation  eleven  included  a  ±  25% 
variation  in  sea-level  pressure  representing  the  sunward  bulge  in  the  atmosphere  and  a  random  variation 


with  a  Gaussian  distribution  of  ±  5%  standard  deviation.  The  orbital  elements  were  chosen  to  represent  a 


realistic  catalog,  with  84%  of  the  catalog  above  800  kilometers  in  altitude. 


Figure  6-28  Simulation  Eleven  Results 

Comparing  this  to  Figure  6-26  shows  that  the  variation  in  atmospheric  density  reduces  the 
accuracy  of  the  estimator.  Although  the  ^nsity  profile  used  in  this  simulation  could  be  made  more 
accurate,  the  important  point  is  that  the  profile  of  the  truth  model  is  far  more  accurate  than  (and  different 
from)  the  estimator  model.  Changing  the  truth  model  to  make  it  more  realistic  would  not  change  the 
results  of  the  simulation.  This  figure  attests  to  the  fidelity  with  which  the  estimator  can  successfiilly 
determine  orbits.  During  the  simulation,  the  average  error  was  as  low  as  two  to  three  kilometers  at  the 
end  of  a  propagate-estimate-update  cycle.  AFSPC  currently  reports  an  average  enor  of  12  kilometers  as 
soon  as  an  estimate  is  made.  This  indicates  the  estimator  developed  in  this  thesis  is  far  more  accurate 
than  that  currently  operational. 

Although  Figure  6-28  indicates  an  improvement  in  estimation  accuracy,  the  results  were  further 
analyzed  to  try  to  find  a  way  to  improve  the  accuracy  even  more.  The  individual  errors  of  all  160 
satellites  were  plotted  as  a  function  of  time.  Figure  6-29  shows  the  results  for  the  first  15  satellites. 
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Figure  6-29  Individual  Satellite  Estimation  Errors 
A  few  low-altitude  satellites  contributed  the  vast  majority  of  the  error  at  each  point  in  the 
simulation.  Satellites  nine  and  fourteen  are  two  of  these  low  altitude  satellites,  with  altitudes  of  perigee  of 
194  and  265  kilometers,  respectively.  The  large  errors  of  these  two  satellites  are  obvious  in  this  figure. 
The  higher  altitude  satellites  were  estimated  with  much  better  accuracy.  The  accuracy  of  satellite  nine 
only  improved  once  during  the  ten  cycles  shown.  This  occurred  when  the  orbit  was  estimated  using  the 
least  squares  batch  estimator  instead  of  the  sequential  Bayes  filter.  These  results  indicate  that  the  low 
altitude  satellites  are  too  difficult  to  estimate  because  of  the  differences  between  the  dynamics  models  of 
the  estimator  and  the  truth  model. 

To  address  this  problem,  a  simulation  was  conducted  using  only  the  five  satellites  with  the  lowest 
altitudes  of  perigee  and  all  nine  tracking  stations.  This  provided  frequent  tracks  of  data  for  the  estimator. 
The  goal  was  to  see  if  a  higher  rate  of  tracking  would  allow  for  more  accurate  orbit  estimation  for  the  low- 
altitude  satellites.  This  simulation  was  conducted  for  1.4  days  and  the  errors  at  each  point  were  recorded. 
Figure  6-30  shows  these  errors. 
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Figure  6-30  Low-Altitude  Satellite  Errors 

This  figure  shows  that,  even  if  the  entire  tracking  network  were  dedicated  to  tracking  only  these 
five  satellites,  the  errors  would  still  average  as  high  as  a  kilometer  or  more.  While  this  seems  to  indicate 
that  further  adjustments  to  the  tracking  priority  would  be  futile,  another  useful  &ct  can  be  taken  fixim  the 
results.  The  errors  of  three  of  the  five  satellites  drof^ied  dramatically  after  the  1.2  day  cycle.  This 
occurred  because  these  three  satellites’  orbits  were  estimated  using  the  least  squares  batch  method  instead 
of  the  Bayes  filter.  This  indicates  that  the  batch  method  is  better  suited  to  the  low  altitude  satellites, 
where  the  changing  atmo^heric  model  makes  the  older  data  much  less  reliable. 

This  analysis  was  used  to  adjust  the  estimator  algorithm  in  two  ways.  First  the  prioritization 
scheme  was  changed  to  place  a  higher  priority  on  satellites  with  tow  altitudes  of  perigee.  Second, 
satellites  which  have  perigee  altitudes  below  350  kilometers  are  flagged  as  low-altitude  satellites  and  their 
oibits  are  only  estimated  using  the  least  squares  batch  method.  These  two  modifications  were  aimed  at 
improving  the  total  accuracy  of  the  estimated  catalog. 

After  these  modifications  were  implemented,  another  simulation  was  conducted  to  show  the 
effects.  All  other  parameters  of  the  simulation  were  identical  to  simulation  eleven.  Figure  6-3 1  shows  the 
results. 
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Figure  6-3 1  Simulation  Twelve  Results 

This  figure  indicates  that  the  modifications  were  successful  in  improving  the  accuracy  of  the 
catalog.  Additionally,  another  set  of  parameters  were  tracked  for  this  simulation.  The  foUowing  figure 
shows  the  average  error  and  the  median  error  in  the  catalog  at  the  end  of  each  cycle.  The  percentage  of 
the  catalog  accurately  estimated  within  one  kilometer  is  plotted,  as  in  the  previous  simulations.  In  this 
plot,  however,  the  percentage  within  100  meters,  within  ten  meters  and  within  one  meter  are  also  shown. 


Figure  6-32  Simulation  Twelve  Error  Trends 
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Figure  6-32  shows  the  trend  in  catalog  accuracy  during  the  simulation.  By  the  end  of  the  four 
day  simulation,  half  of  the  satellites’  positions  were  predicted  within  63  meters.  This  represents  a  very 
significant  improvement  over  the  current  method.  Because  the  results  were  better  with  the  increased 
likelihood  of  batch  estimations  for  low  altitude  satellites  (simidation  twelve)  over  the  previous  result 
(simulation  eleven),  another  simulation  was  conducted  to  see  if  increasing  the  tendency  toward  batch 
estimations  for  all  satellites  would  increase  the  accuracy  further.  To  investigate  this,  the  algorithm  was 
modified  so  that  the  batch  method  was  implemented  after  any  unsuccessful  sequential  estimation  attempt 
with  at  least  three  new  tracks  of  data.  Figure  6-33  shows  the  results  of  this  simulation. 


day 

Figure  6-33  Simulation  Thirteen  Results 


The  results  were  better  than  simulation  twelve  but  only  marginally.  Nevertheless,  this  simulation 
represents  the  best  accuracy  obtained  during  this  thesis  with  a  realistic  catalog  and  a  realistic  truth  model. 
Consistently,  more  than  85%  of  the  catalog  was  estimated  within  one  kilometer  and  the  median  error 
averaged  less  than  50  meters. 

The  third  atmospheric  effect  modeled  is  that  of  a  significant  solar  flare.  The  effect  of  this  flare 
on  the  estimation  process  was  analyzed  for  the  final  simulation.  Figure  6-34  shows  the  atmospheric 
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densi^  profile  modeled  for  this  simulation.  It  includes  a  day-night  variation  of  ±  25%  and  a  solar  flare 
variation  of  a  factor  of  two  in  sea-level  pressure. 


altitude  (kilometers) 

Figure  6-34  Atmospheric  Density  Variations 

As  shown  above,  the  d^-night  variation  in  sea-level  pressure  of  ±  25%  results  in  a  change  in 
density  of  ±  1 1%  at  200  kilometers.  At  this  same  altitude,  the  increase  in  sea-level  pressure  Ity  a  factor  of 
two  fixnn  the  solar  flare  results  in  an  increase  in  density  Ity  a  factor  of  two.  To  model  the  flare  in  the  truth 
model  (tynamics,  the  sea-level  pressure  was  instantaneously  doubled,  and  then  exponentially  returned  to 
normal  during  the  following  day.  This  approximately  models  the  density  changes  at  altitude  resulting 
from  a  large  solar  flare. 

Simulation  fourteen  shows  what  would  happen  to  this  estimation  process  if  a  solar  flare  of 
significant  magnitude  were  to  affect  an  increase  in  atmospheric  density.  This  scenario  is  quite  plausible 
and  rei»esents  a  significant  obstacle  for  an  estimation  algorithm  to  overcome.  A  large  solar  flare  will 
dramatically  and  unpredictably  increase  the  atmospheric  density  at  a  particular  altitude.  This  causes  a 
major  increase  in  atmospheric  drag  on  low-altitude  satellites  causing  large  orbital  perturbations  and  even 
premature  reentry.  The  solar  flare  began  at  the  end  of  the  second  day,  and  returned  to  normal  near  the 
end  of  the  third  day.  This  simulation  was  expected  to  demonstrate  the  significant  impact  which  a  solar 
flare  wiU  have  on  the  estimation  accuracy.  Figure  6-35  shows  the  results  of  this  simulation. 
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Figure  6-35  Simulation  Fourteen  Results 

As  expected,  this  figure  shows  that  the  solar  flare  dramatically  impairs  the  estimator’s  ability  to  accurately 
predict  orbits.  Even  in  this  worst  case  scenario,  the  track  compression/G/ofto/  Estimate  process  still 
recovers  to  outperform  the  current  method. 


6.8  Summary 

In  summary,  the  simulations  showed  that  the  track  compression/Global  Estimate  method  of  orbit 
estimation  will  work  on  a  full  size  catalog  and  that  the  results  are  much  better  than  current  practice.  The 
results  showed  that  different  types  of  orbits  produce  a  wide  variation  in  acciuacies,  with  low  altitude  orbits 
being  the  most  difficult  to  estimate.  This  is  both  because  the  low  altitude  <fynamics  are  stronger,  meaning 
that  any  error  will  grow  quickly,  and  because  the  low  altitude  satellites  are  subjected  to  the  impredictable 
atmosphere,  making  an  estimator’s  dynamics  model  less  valid. 

The  results  also  showed  that  the  covariance  matrices  are,  in  general,  overly  optimistic  in  the 
prediction  of  the  amount  of  error  in  any  estimate.  Several  attempts  were  made  at  making  the  covariance 

I 

matrices  a  better  indicator  of  estimator  performance.  The  application  of  fading  memory  coefficients  did 
not  significantly  improve  the  results.  One  simulation  which  did  improve  the  results  placed  the  Global 
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Estimate  dynamics  model  and  truth  model  in  perfect  agreement.  When  the  more  realistic  atmospheric 
effects  were  added,  the  covariance  accuracies  were  again  optimistic.  This  is  one  area  which  requires 
further  evaluation.  A  more  accurate  fading  memory  matrix  with  unique  values  for  all  elements  should  be 
developed.  The  best  values  potentially  are  altitude-specific. 

Finally,  phase  three  demonstrated  that  the  track  compression/G/of)u/  Estimate  method  would 
produce  accuracies  measured  in  meters  or  tens  of  meters  rather  than  kilometers.  The  method  requires 
modest  computer  time.  Time  trials  show  that  a  133  MHz  personal  computer  ruiming  in  a  Windows 
environment  could  handle  the  real-time  processing  of  approximately  320  to  400  satellites  using  the 
algorithm  developed  in  this  thesis. 
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VII.  Conclusions  and  Recommendations 


7J  Summary 

This  thesis  presents  two  successful  methods  of  compressing  tracks  of  observation  data  into  state 
vectors  and  covariance  matrices.  It  also  presents  an  estimator  which  determines  orbits  using  track 
compression  results  as  input  data  and  shows  the  applicability  of  this  process  to  the  satellite  tracking 
mission  by  simulating  its  use  on  large-scale  catalogs. 

7.2  Conclusions 

To  conclude,  the  track  compression  methods  were  successful.  The  Taylor  series  method 
successfully  compressed  every  arc  to  which  it  was  applied  Because  of  the  requirement  to  split  each  track 
into  multiple  arcs,  this  method  was  abandoned  in  favor  of  the  integrator  method.  The  integrator  method 
successfully  compressed  over  one  hundred  thousand  tracks  of  data  taken  from  thousands  of  different 
orbits.  This  represents  a  success  rate  of  greater  than  99.99%.  Additionally,  the  method  showed  that  the 
orbits  could  be  determined  even  though  the  track  compression  dynamics  model  was  much  less  accurate 
than  the  truth  model.  In  simulation  ten,  the  track  compression  (fynamics  model  did  not  include  air  drag, 
yet  the  orbit  estimator  was  still  able  to  use  the  track  compression  results  to  predict  orbit  positions  within 
one  kilometer  for  more  than  99%  of  the  catalog.  This  shows  that  the  state  vectors  and  covariance  matrices 
retain  enough  of  the  accuracy  of  the  original  arcs  that  a  sufficiently  accurate  orbit  estimator  can  correctly 
determine  orbits  with  this  greatly  reduced  amount  of  data. 

Additionally,  the  simulations  showed  that  the  Global  Estimate  algorithm  would  produce  a 
significant  improvement  in  catalog  accuracy  over  current  methods.  Today,  the  SCC  reports  orbit 
estimations  with  a  one  sigma  error  of  twelve  kilometers.  This  indicates  that  about  twenty  percent  of  the 
orbit  determinations  are  within  twelve  kilometers.  The  simulation  results  indicate  that  the  methods 
developed  in  this  thesis  would  result  in  twenty  percent  of  the  orbit  predictions  within  ten  to  twenty  meters. 
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This  represents  an  improvement  of  two  to  three  orders  of  magnitude  over  the  orbit  estimation  models 
currently  in  use  the  SCC. 

There  are  also  two  reasons  which  indicate  that  the  actual  accuracies  may  be  much  better  than 
those  predicted  hy  the  simulations.  First,  the  latest  simulations  were  conducted  on  a  catalog  with  84%  of 
the  semi-major  axes  above  800  kilometers.  A  more  accurate  catalog  would  probably  have  84%  of  the 
altitudes  of  perigee  above  800  kilometers.  [4]  The  estimation  process  conducted  in  the  latest  simulations 
focused  much  attention  on  trying  to  estimate  orbits  which  were  low  enough  to  have  significant  drag 
effects.  These  satellites  produce  the  majority  of  the  error  of  the  catalog.  With  the  actual  catalog,  most 
satellites  are  not  nearly  as  influenced  hy  air  drag  as  those  modeled  in  this  thesis.  Therefore,  the  estimator 
would  spend  less  time  trying  to  estimate  the  low  altitude,  difBcult  satellites  and  would  produce  more 
accurate  predictions  of  the  higher  altitude  satellites. 

Secondly,  the  one  sigma  accuracy  reported  by  the  SCC  applies  at  the  instant  that  a  satellite’s 
orbit  is  determined.  This  process  is  called  filtering.  [7]  The  one  sigma  accuracy  reported  during  the 
simulations  is  taken  at  an  epoch  corresponding  to  the  end  time  of  each  cycle.  This  could  be  several  orbits 
after  the  reference  time  of  the  orbit  determination.  This  process  is  called  prediction.  [7]  State  prediction 
is  much  less  accurate  in  orbital  mechanics  than  filtering  because  even  the  best  estimate  of  an  orbit  will 
lose  accuracy  with  time.  Also,  if  a  limited  number  of  observations  are  available  for  a  satellite,  the  opoch 
of  the  prediction  could  be  far  from  the  time  of  the  observations,  meaning  that  the  estimator  is  trying  to 
predict  where  the  satellite  will  be  on  one  side  of  its  orbit  based  on  information  from  only  the  opposite  side. 
Both  of  these  reasons  point  to  an  improvement  potentially  greater  than  the  two  to  three  orders  of 
magnitude  stated  above. 

7.3  Recommendations 

This  thesis  showed  that  the  use  of  the  track  compression/G/ofea/  Estimate  method  of  orbit 
estimation  for  the  operational  satellite  tracking  mission  should  be  considered.  The  computational 
requirements  of  this  new  method  are  modest  enough  to  be  implemented  by  the  SCC.  It  would  greatly 
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improve  the  acxniracy  of  orbit  determinations  which  would  result  in  a  reduced  routine  radar  tracking 
workload,  fewer  collision-avoidance  maneuvers  for  maimed  space  flights,  more  accurate  over-flight 
predictions,  faster  acquisition  for  ground-based  observation  or  targeting  of  satellites  and  smaller 
rendezvous  fuel  budgets. 

There  are  several  areas  which  remain  to  be  investigated.  These  deal  with  the  operational 
implementation  of  the  algorithm.  First,  the  long  term  behavior  of  the  estimator  needs  to  be  evaluated. 

This  will  require  extended  computer  simulation  time.  A  Monte-Carlo  analysis  should  be  conducted  to 
determine  appropriate  values  for  each  element  of  the  fading  memory  coefficient  matrix.  The  optimal 
values  for  these  elements  will  likely  be  functions  of  altitude. 

A  number  of  different  techniques  were  employed  to  increase  the  reliability  of  the  code.  Each  of 
these  should  be  investigated  to  determine  any  adverse  effects  caused  by  their  implementation.  These 
include  the  exclusion  from  the  estimation  process  tracks  of  data  which  have  less  than  40  observations  and 
the  exclusion  of  tracks  with  less  than  240  observations  from  use  as  the  reference  trajectory  during  the 
estimation  process.  The  ballistic  coefiBcient  was  carefully  monitored  during  the  convergence  process.  The 
algorithm  had  to  give  a  reasonable  ballistic  coefBcient  at  each  iteration,  or  else  it  would  be  stopped 
Additionally,  the  algorithm  had  to  converge  on  a  tighter  set  of  ballistic  coefficient  criteria.  Finally,  a 
check  was  made  at  each  iteration  to  see  if  this  parameter  was  diverging. 

During  the  propagation  of  the  state  vector  at  each  iteration  of  the  estimation  process,  the  altitude 
and  energy  of  the  oibit  were  checked  to  see  if  the  orbit  was  being  propagated  too  low  into  the  atmosphere, 
or  if  the  orbit  had  escaped  The  former  would  be  an  indication  that  the  satellite  was  very  near  reentry 
while  the  later  is  a  sign  of  divergence.  Additional  divergence  checks  were  also  used. 

The  track  compression  algorithm  would  give  up  if  it  failed  to  converge  within  15  iterations.  The 
same  limit  was  placed  on  the  Global  Estimate  algorithm.  The  residuals  were  checked  after  convergence 
to  see  if  the  result  was  reasonable.  If  the  algorithm  detected  that  a  covariance  matrix  diagonal  element 
had  become  negative,  it  was  assumed  that  this  was  the  result  of  a  numerical  precision  limitation  and  that 
estimate  was  postponed  until  more  data  was  available.  Also,  if  an  oibit  could  not  be  successfully 
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estimated  using  the  Bayes  sequential  method,  and  30  or  more  tracks  of  new  data  were  available,  the  new 
data  would  be  fed  into  a  least  squares  batch  estimator.  For  the  last  two  simulations,  this  was  changed  to 
three  or  more  new  tracks  of  data.  If  the  satellite  were  flagged  as  a  low  altitude  satellite,  only  batch 
updates  were  performed. 

Finally,  the  satellites  were  prioritized  in  part  based  on  the  size  of  the  covariance  matrices,  but 
also  on  the  nmnber  of  successive  &iled  estimation  attempts  and  the  altitude  of  the  satellite.  All  of  these 
methods  of  modifying  the  basic  theoretical  algorithm  produced  either  more  robust  code,  more  accurate 
results  or  both.  The  optimal  means  of  implementing  each  of  these  modifications  should  be  investigated 
both  to  determine  if  there  is  a  more  effective  way  of  modifying  the  algorithm  or  if  the  change  adversely 
affects  the  process.  Some  of  these  may  require  more  processing  time  than  they  are  worth.  The  multitude 
of  combinations  of  the  above  adjustments  makes  this  optimization  stufy  far  too  time  consuming  to  be 
attempted  in  this  thesis. 

Additional  variables  which  were  not  investigated  significantly  in  this  thesis  include  varying  the 
length  of  time  for  each  cycle  and  estimating  orbits  using  altitude-specific  estimation  algorithms.  Also,  a 
method  of  prioritizing  satellites  which  looks  at  which  parts  of  the  orbit  have  not  been  observed  could  be 
considered.  Finally,  the  covariance  matrix  for  each  satellite  could  be  propagated  along  with  the  state 
vectors  during  the  observation  portion  of  the  cycle,  and  reprioritized  at  each  time  step.  This  may  be  a 
more  optimal  prioritization  scheme.  Because  the  change  in  the  covariance  matrix  due  to  new  data  is  only 
a  function  of  the  accuracy  of  the  new  data,  the  covariance  matrix  could  be  altered  after  a  new  track  is 
compressed  even  before  an  estimate  is  performed  These  and  many  other  methods  of  optimizing  the 
fuioritization  algorithm  could  be  explored 

By  exploring  the  effects  of  the  robustness  modifications,  examining  potential  improvements  not 
investigated  in  this  thesis  and  optimizing  the  prioritization  algorithm,  the  track  compression/G/oAu/ 
Estimate  method  of  orbit  estimation  can  be  improved  and  prepared  for  operational  use. 
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The  results  indicate  that  a  significant  imiHovement  in  estimation  accuracj’  would  be  seen  if 
AFSPC  were  to  move  to  this  method.  The  use  of  this  method  for  the  operational  satellite  tracking  mission 
should  be  considered. 
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Appendix  A  Atmospheric  Model 


For  the  truth  model,  Global  Estimate  and  Simulate  portions  of  this  thesis,  an  atmospheric  model 
was  required  The  model  provided  the  density  and  density  gradient  as  a  function  of  altitude.  The 
atmospheric  model  used  is  taken  from  Regan.  [15]  It  is  derived  from  a  combination  of  two  standard 
atmospheric  models.  These  models  essentially  provide  a  description  of  the  variation  of  temperature  with 
altitude  and  the  relationship  between  density  and  temperature.  The  models  divide  the  atmosphere  into 
altitude  regions  called  strata.  The  relationship  between  temperature  and  altitude  is  different  for  each 
strata. 

The  first  model,  U.S.  Standard  Atmosphere  1976,  provides  seven  strata  between  sea-level  and  86 
km  of  altitude.  Within  each  of  these  strata,  the  temperature  either  varies  linearly  with  altitude  or  the 
strata  is  isothermal.  Above  86  km,  four  more  strata  serve  to  define  the  atmosphere  up  to  an  altitude  of 
1000  km.  One  of  these  regions  defines  an  elliptic  variation  of  temperature  with  altitude.  A  second  region 
defines  an  exponential  variation.  To  keep  the  model  simpler,  Regan  recommends  that,  above  86  km  a 
1962  model  be  used.  This  model  defines  thirteen  regions  between  86  and  700  km,  all  with  linear 
temperature-altitude  relationships.  Attempts  were  made  to  prevent  the  propagation  of  a  satellite  into 
regions  of  the  atmosphere  defined  by  the  1976  model  but  the  code  describing  these  strata  was  included  for 
the  occasion  when  those  altitudes  had  to  be  evaluated 


The  density,  p,  for  a  given  altitude  is  determined  from  the  following  two  equations 


where  i  is  the  region 

Pi  is  the  density  at  the  base  of  the  region 
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go  is  the  sea-level  gravitational  acceleration 

Z  is  the  altitude  in  kilometers 

Zi  is  the  altitude  at  the  base  of  the  region 

R  is  the  atmospheric  gas  constant,  given  1^  R  =  R*/Mo 

R*  is  the  imiversal  gas  constant 

Mo  is  the  sea-level  molecular  weight  of  air 

Tmi  is  the  molecular  temperature  at  the  base  of  the  region 

b  =  2/RE 

Re  is  the  radius  of  the  Earth 

Lzi  is  the  thermal  lapse  rate,  the  linear  constant  of  variation  of  temperature  for  the  region 
The  following  table  provides  the  values  used  for  many  of  these  constants. 

Table  A-1  Atmospheric  Constants  [15]  , 


1  ^ 

R* 

8.31432  X  10^  (kg-mole)  ‘ 

1  ^ 

28.964  kg/(kg-mole) 

1 

6.3781  X  10®  m 

The  values  for  the  different  strata  are  provided  in  the  following  two  tables. 


Table  A-2  1976  Standard  Atmosphere  [15] 


layer  index 

geometric  altitude 
km 

molecular  temperature 

K 

lapse  rate 
K/km 

0 

0.0 

288.15 

-6.5 

1 

11.0102 

216.65 

0.0 

2 

20.0631 

216.65 

1.0 

3 

32.1619 

228.65 

2.8 

4 

47.3501 

270.65 

0.0 

5 

51.4125 

270.65 

-2.8 

6 

71.8020 

214.65 

-2.0 

7 

86.0 

186.946 
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Table  A-3  1962  Standard  Atmosphere  [15] 


layer  index 

geometric  altitude 
km 

molecular  temperature 

K 

lapse  rate 
K/km 

7 

86.0 

186.946 

1.6481 

8 

100.0 

210.65 

5.0 

9 

110.0 

260.65 

10.0 

10 

120.0 

360.65 

20.0 

11 

150.0 

960.65 

15.0 

12 

160.0 

1110.65 

10.0 

13 

170.0 

1210.65 

7.0 

14 

190.0 

1350.65 

5.0 

15 

230.0 

1550.65 

4.0 

16 

300.0 

1830.65 

3.3 

17 

400.0 

2160.65 

2.6 

18 

500.0 

2420.65 

1.7 

19 

600.0 

2590.65 

1.1 

20 

700.0 

2700.65 

These  two  tables,  along  with  equations  ( A 1 )  and  ( A.2 )  are  used  to  determine  the  density  for  a  given 
altitu(te.  Equation  ( A.  1 )  applies  to  the  isothermal  regions  (where  Lzi  =  0)  while  equation  ( A.2 )  ai^lies 
to  the  linear  regions. 

To  determine  the  density  gradient,  (the  derivative  of  the  density  with  respect  to  altitude), 

equations  ( A 1 )  and  (  A.2 )  ate  differentiated.  Recognizing  that 

9Z  dZ 
3r  3alt 

leaves 
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0*1  C2l 


(A.3) 


dp  pjgo 

dr  RT^ 


[b(z-Zj)“l]expj 


go(Z-Zi)' 

RTm, 


.i-|(z-Zi)] 


for  the  isothermal  regions  and 


L. 

J 

f  So  1 

RLzi  { 

'Zk.„z  )1 

[ULziJ 

- +  1+  D| 

go  \ 

gpb 

RLa 


^  go 


RL 


ZiJ 


RL, 

So 


+  l  +  b 


f  T 

^Mi  -7 

T 

V^zi 


A 


for  the  linear  regions. 


(A.4) 
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Appendix  B  Geopotential  and  the  Equations  of  Motion 


The  purpose  of  this  appendix  is  to  present  the  oblate  Earth  model  of  the  geopotential  and  to 
develop  equations  of  motion  for  a  satellite  imder  the  influence  of  this  potential.  These  equations  of  motion 
were  used  in  the  truth  model  to  propagate  the  satellites’  orbits.  Since  the  accelerations  due  to  the  Earth’s 
geopotential  are  only  functions  of  the  satellite’s  position,  the  following  simplified  notation  will  be  used 

R  =  [x  y  zf 

and  the  equations  of  motion  will  take  the  form 

’d'x' 

dt" 

d^R  d^y 
dt^  "  dt^ 
d^z 
dt' 

The  geopotential  model  presented  here  does  not  include  sectoral  or  tesseral  harmonics. 

The  aspherical  potential  of  the  Earth,  modeling  zonal  harmonics  through  Je,  is  provided  1^ 
Escobal  as  [5] 

«  =  —  l+^^(l-3sin’S) 

H — ^^^{3-5sin^5)sin5 

R  4  j 

^-r-^(3-30sin^  S  +  SSsin^S)  (B.l) 

8r 

TOsin' 5  +  63sin^  8)sin5 

R  6  j 

+ — ^-^(5  - 1 05  sin' 6  +  3 1 5  sin^  5  -  23 1  sin®  5) 

16r®  ' 
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where  G  is  the  gravitational  constant 
m  is  the  mass  of  the  Earth 

is  the  equatorial  radius  of  the  Earth 

r  is  the  distance  from  the  Earth’s  center 
sin  8  =  z/t 

and  the  zonal  harmonic  coefficients  are  taken  from  the  following  table  [5] 


Table  B-1  Coefficients  of  Gravitational  Harmonics 


h 

+ 

1082.28 

±0.3x10^ 

h 

- 

2.3 

±0.2x10"® 

Ja 

- 

2.12 

±0.5x  10"® 

h 

- 

0.2 

±0.1x10"® 

h 

1.0 

±0.8x  10"® 

Using  canonical  units  allows  the  substitution  s  Gm  =  1  and  =  1.  The  acceleration  due  to  gravity 
about  a  planet  is  given  as  the  gradient  of  the  potential,  or 


d^R 

dt' 


=  Vd> 


(B.2) 


and  these  terms  are  given  as 


dx 


+■^^(3— 7sin*  8)sin5 

Y 

-  ^^(3  -  42  sin"  5  +  63  sin''  5) 
o  r 

3  J 

-  --^(35  -  210  sin"  5  +  23 1  sin''  5)  sin  5 

+ ^^(35  -  945sin"  5  +  3465  sin''  5  -  3003  sin" 


(B.3) 


dy  dx  J 
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dz 


“7  ^  + 

+  ^^(6-7sin^  5)  sin  5 

- (1 5  -  70  sin'  5  +  63  sin^  5) 

8  r 

- ^  {1 05  -  3 1 5  sin'  5  +  23 1  sin'*  5)  sin  5 

+  —^(245  -  2205  sin'  5  +  4851  sin'*  5  -  3003  sin"  5 
16  r"  ^ 


■^r'Ur^  8r"j 


Combining  these  yields  the  final  equations  of  motion  for  a  satellite  about  an  oblate  planet. 


(B.4) 
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Appendix  C  Taylor  Series  Track  Compression  State  Transition  Matrix 


This  appendix  presents  the  state  transition  matrix,  O ,  used  in  the  Taylor  series  method  of  track 
compression.  The  state  transition  matrix  relates  a  change  in  the  elements  of  the  position  vector  at  time,  t, 
to  a  change  in  the  reference  state,  X ,  at  time,  to. 


aR(t) 

9X(0) 


(C.l) 


where  the  following  notation  is  used 


R(t)  =  R  = 


R,(t) 

R2(t) 

.RsO). 


X(0)  = 


R, 

R2 

R3 

V, 

% 

.V3. 


Using  this  notation,  the  state  transition  matrix  becomes. 


,  aR(t) 

dX(0) 


3R.(l) 

ax,(o) 

3R.(t) 

ax,(o) 

3R,(t) 

ax,(o) 


aRi(t) 

3X,(0) 

aRi(t) 

ax,(o) 

aRjW 

ax,(o) 


aR,(t)  aR,(t)  aR,(t)  aR,(t) 


ax,(o) 

ax.(o) 

ax,(o) 

ax,(o) 

aR,(t) 

aR,(t) 

aR,(t) 

aR,(t) 

ax,(o) 

ax,(o) 

ax,(o) 

ax,(o) 

3R,(t) 

aR,(t) 

aR,(t) 

aR,(t) 

ax,(o) 

ax,(o) 

ax,(o) 

ax,(o) 

To  determine  these  partial  derivatives,  R(t)  as  a  function  of  X  must  first  be  written.  This  function  is 
developed  in  section  3A.2.2  as  equation  ( 3.15 )  and  is  reprinted  on  the  following  page. 
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R 


Ri  ^  3  R, 

R,+V,At-— VAt'-- 

I  1  2P  4r* 


1-5- 


R 


2^ 


R2  +  V2  At  ‘ 


I2. 

2r" 

R3 

R3-l-V3At--^ 

3  3  2r^ 


At' 


3J2R2 


At'- 


4r* 
3J3R3 


2^ 


At'  + 


1 

r  ''1 


6r' 


At' 


y 


1-5- 


R 


4r 


3-5- 


r 

R 


) 
2  ^ 


,  1  fSRjf  ^  3 

At'-I--^  — —-^2  At' 

6r'l  r  M 


At'  +■ 


1  r3R3f 


6r' 


■V3  At' 


(3.15) 


where  r  =  ■y^Rj'~+R^'”+R^ 

dt  r 

At  =t-to 


J2  =  1082.28  X  10-* 

Replacing  the  individual  terms  of  equation  (3.15)  with  another  shorthand  notation  will  allow  the 
components  of  the  partial  derivatives  to  be  developed  piecemeal.  The  equation  is  reprinted  below  in  the 
same  form  to  allow  an  easy  understanding  of  which  notation  corresponds  to  each  term. 


R  = 


Ajj  -f-  At  +  C^At'  -l-  Djf  At'  -l-  Ex  At' 
Ay  +  By  At  +  CyAt'  +  Dy  At'  +  Ey  At' 
^  A2  4-  B^  At  4-  C^A  t'  4-  At'  4-E2At' 


(C.2) 


As  an  example,  the  first  partial  derivative  of  equation  (3.16)  becomes 

aR,(t)  BAx  aBx  acx  ,  aox  ,  aEx 
ax,(o)  ax,(o)  ax,(o)  ^  ax,(o)  ax,(o)  ax,(o)  ^  ' 


and  is  abbreviated 

aR,(t)  ,  ,  3 

a X  (0)  ~  ^ ®x  1 A 1 4- Cx  1 A t  4-DxlAt  4-Ex  1  At  (C.4) 

Therefore,  determining  the  eighteen  components  of  equation  (3.16)  becomes  a  process  of 
determining  the  five  partial  derivatives  of  equation  ( C.3  )  for  each  of  the  eighteen  equations  which  have  a 
form  identical  to  ( C.3  ).  These  90  terms  are  denoted  1  through  E^  6 .  These  will  be  much  simpler 
to  obtain  and  describe  than  if  all  five  terms  were  evaluated  as  one.  In  many  cases,  they  are  trivial. 
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Use  of  a  simplified  notation  will  also  be  seen  in  the  extensive  application  of  the  chain  rule.  The 
partial  derivatives  of  certain  terms  will  occur  frequently  and  will  not  be  evaluated  each  time.  Instead,  the 
term  will  be  indicated  in  the  appropriate  equation  and  evaluated  below.  In  the  actual  computer  code,  this 
process  occurs  in  the  reverse  order.  The  terms  are  evaluated  and  stored  beforehand,  and  then  retrieved 
frequently  during  the  evaluation  of  the  appropriate  equations. 

The  remainder  of  this  aj^ndix  is  devoted  to  presenting  each  of  the  90  partial  derivatives. 


Axl  =  l 

Ax2  =  0 

0 

it 

cn 

X 

< 

Ax4  =  0 

Ax5  =  0 

Ax6  =  0 

0 

It 

>< 

< 

Ay2  =  1 

Ay3  =  0 

Ay4  =  0 

Ay5  =  0 

Ay6  =  0 

A,1  =  0 

Az2  =  0 

Az3  =  l 

Az4  =  0 

Az5  =  0 

Az6  =  0 

Bxl  =  0 

Bx2  =  0 

Bx3=0 

Bx4  =  l 

Bx5  =  0 

Bx6  =  0 

By1  =  0 

By2  =  0 

By3  =  0 

By4  =  0 

By5  =  1 

By6  =  0 

B^1  =  0 

B,2  =  0 

B23  =  0 

Bz4  =  0 

Bz5  =  0 

Bz6  =  l 

r  1-  -1 

,  3R1  ar 

r  2  3R,  ar 

r  T  3R,  ar 

^  2r^ 

2  aR, 

^^^“2  aR, 

aR, 

Cx‘»  =  0 


Cj5  =  0 


Cx6  =  0 


C  1  =  1^-^ 
’  2  r*  3R, 


C  2  =  ^+1^-^ 

2r=  2  r'  aR, 


3R,  ar 


r  3=- 

2  dR, 


Cy4  =  0 


Cy5  =  0 


Cy6  =  0 


C,l  = 


3^_8r_ 

2  r'*  BR, 


Cz2  = 


3^_£L 

2  dR, 


-1  3R3  Br 

Cz3  =  :rT+- 
2r 


2  r'*  BR, 


Cz4  =  0 


Cz5=0 


C^6  =  0 


r 


Dxl  = 


2  ^ 


Ri-5Ri 


I5J2 


4  r‘  aR, 


R, 


R3^  ar 


1-5^  +  lOR,^^ 


Dx2  = 


R,-5R, 


R3M15J2  3J2r,.„  R3' 


4  r®  aR,  4  r’ 


ill. 

4r* 


lOR, 


aR 
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f  Rs'llsJz  3J2 

Dx3  =  [R.-5R,-^Jf7i-^-47 


Dx4  =  0 


y  ■  •  •'^'3  r 

Dx5=0 


R3'  1 

Dx6  =  0 


Dy2  = 


Dy1  = 

r 


_cp  R3'1i5J2  3J2 

2  ^^2  J.2  J  4  J.6  4  J.J 

R3^^15J2  3r 

p  ^  _L£.— £ _ _ 

«-2  5K3  ^2  1  4 


2 1  R3' 

10R2  J.3 


R,^  R3"  dr 

1"5  J.2  +IOR2  ^3 


2; 


\lb. 

4r^ 


Dy3  = 
Dy4  =  0 

Dzl  = 


R/] 
^^2  2 
^  J 


il£2._^_3l2.f,f.p  R3'  ^ 

4r'^aR3  r^  aR3j 


Dy5  =  0 


r^  dRsj 

Dy6  =  0 


R30 

^  -n  r  i 

15J2  ar  15J2 

fR,’  Sr  7 

^^3  ^  2 

V  ^  ) 

4  r"  aRi  2  r’ 

lr>  SrJ 

D  2-f3R 

U^z-^JK,  :>  ^2  J  4  J.6  2  r' 


^R3^  ar  ^ 


p  aR. 


■2j 


Dz3  = 

V 

Dz4  =  0 

_i  ar 
Evi  =  .  i 


R3M15J2  dr  _R3^  _R3^  arl3J2 

3R3-5-J-|f7^-[3-l5-7^^l^— ^J47 


Dz5  =  0 


2  =  2?' r  "  ir* 


-1  ar 

'  2r"  aR3 


r  ^ 


(  r 

3-+3 


( 


R, 

dr 

r 

■^aR, 

Ri 

dr 

I 

r 

aR2 

Ri 

af 

1  -  -  - 

r 

aR3 

■aR3j4r 

Dz6  =  0 

f  ar  ^ 

aR,, 

f  ar  ^ 

f  ar 
‘r^  aR3; 


1  R,  af  ^ 
Ex4  =  7V  3— ttt-I 
^  6r^V  r  av,  ^ 
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E  5  =  —^ 
^  6r^ 


Ev6  = 


r  R.  3^ 

J  T  r  -I 


6f 


Evl  = 


-1 
iv^  aR 

1  ar 


r  av. 


2  y 
^ 


R,  af 
X  dW..  ^ 


3  J 


f 


Ev2  = 


Ev3  = 


2r"  aR^  V  ^  r 
-1  ar 


3R,--V,  + 


6r' 


^"■2r^aR3 


1  LR2 


6r^ 

1’  r  3V, 

1  1 

'  R,  3f 

6rM 

y  r  av. 

1 

f  R,  3f 

6r" 

r  f  3V3 

-1  ar 

"''■2r^  aR3 


E,3=: 


[3R3i-V3)+^ 


R2 

af 

► 

r 

aR, 

af 

1  - . 

r 

^aR^ 

,R2 

dr 

1 

r 

^aR3 

■> 

--1 

J 

\ 

-1 

y 

-1 

J 

R3 

dr 

1 

r  aR, 

R3 

dr 

r  "^aR^ 

R3 

dr 

r 

^dR, 

r  ar 


^^2  2 


r  aR, 


,  ±_?I  "l 

r'  aR 


2J 


-3R 


f  ar 

2  r'  aR3 


-3R 


r  ar 
3  -2  aR 


^ 


-3R 


r  w  i'-i  y 

f  ar  ^ 

'  aR, 


-3R 


f  ar 


^ 


aR 


3  y 


1  r  R3  af  " 

£74  =  -^  3— T— -1 

^  6rn  r  av,  ^ 


E  5  =  ^^ 
^  6r' 


R3  af  ^ 

‘  3Tr-i 

r  aVj  j 
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Ez6  = 


6r' 


R,  8f 
r  BY, 


where  r  = 


ar  R, 

ar  R2 

Br 

aR,  “  r 

aR2  “  r 

aR, 

r 

.  d  RjV, +R2V2+R3V3 

r 

Br  V,  fR, 

af  V2  fR2 

af 

V3 

aR,  ■  r  ”  r' 

aR,  -  r  “  r' 

aR, 

”  r 

af  R, 
av,  ”  r 

II 

dr 

av, 

r 

Inserting  these  substitutions  into  the  90  partial  derivatives,  inserting  these  90  partial  derivatives  into  the 
eighteen  equations  of  the  form  ( C.4 ),  and  finally  inserting  these  eighteen  equations  into  equation  (3.16) 
yields  the  final  form  of  the  state  transition  matrix,  O ,  for  the  Taylor  series  method  of  track  compression. 
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Appendix  D  Track  Compression  H  Matrix 


This  appendix  develops  the  H  matrix,  described  in  section  3.4.4. 1 .  This  matrix  relates  the 


change  in  the  observation  relation,  G ,  with  a  change  in  the  position  vector,  R . 


3Gj 

aG, 

aG, 

aR, 

aR, 

dR, 

ac, 

ac, 

aR, 

aR, 

aRs 

ac, 

ac, 

aG, 

aR, 

aR, 

aR, 

The  observation  relation,  G ,  is  given  the  following  equations 


G  = 


range(R,t) 
azimuth  (R,t) 


elevation(R,t) 


(3.28) 


(3.17) 


range  —  |  ~  VPs  "^Pe  '^'Pz 


(3.11) 


azimuth  =  Jt  -  tan  ’ 


Pe 

Ps 


elevation  =  tan"'  ,  — 

a/Ps  "^Pe 

Because  the  partial  derivative  of  G  must  be  taken  with  respect  to  R ,  and  G  is  expressed  in  terms  of 
>  the  relations  between  the  various  terms  are  used  to  create  a  chain  of  derivatives  as  follows. 

3G  dG  3p|5Ez  ^Pdk  .p. 

"-aR-8fe  api„,  aR  ' 

where  G  is  given  in  terms  of  fi-py  by  equations  (3.17)  and  (3.11)  above,  p^E2  Pdk  related 
through  equation  ( 3.23  ),  and  ppK  and  R  are  related  through  equation  (3.18),  both  on  the  next  page. 


D-1 


(3.23) 


PsEZ  ~  ®  PuK 


PdK  ~  ^(0  ^she 


The  three  terms  of  equation  ( D.  1 )  are  now  determined  in  reverse  order.  From  equation  (3.18), 


1 


^Pdk 

dR 


=  1  = 


0 


0 


0  0* 
1  0 
0  1 


From  equation  (  3.23  ), 


where 


^PsEZ 

^Pdk 


=  D 


D  = 


cos(LST)  cos(colat) 
-  sin(LST) 
cos(LST)  sin(colat) 


sin(LST)  cos(colat)  -  sin(colat) 
cos(LST)  0 

sin(LST)  sin(colat)  cos(colat) 


The  third  derivative. 


3G 

^PlsEZ 


,  will  now  be  developed  from  equations  (3.17)  and  (3.11). 


3  range 

Ps 

3ps 

+  Pe'+Pz' 

3  range 

Pe 

^Pe 

VPs^ 

+  Pe'+Pz' 

3  range 

Ps 

5ps 

VPs' 

+  Pe'+Pz' 

3  azimuth 

-Pe  Ps“' 

3ps 

1 

ffp/T 

A  /Ps) 

(3.18) 


(D.2) 


(D.3) 


(  3.22  ) 


(D.4) 


D-2 


d  azimuth 
9pE 


1  + 


9  azimuth 


9 elevation  -  pz  ps  (ps^  +  Pe^  )  ^ 

^Ps  "  Ps'+Pe'+Pz' 

9 elevation  -  pz  Pe  (Ps^  +  Pe^  ) 

9Pe  "  Ps'+Pe'+Pz' 


9  elevation  VPs^+Pe^ 

^Pz  (Ps^ +Pe^ +Pz^) 


These  nine  equations  are  inserted  into  the  following  matrix. 

9  range 

9  range 

9  range 

9G 

3ps 

9  azimuth 

3Pe 

9  azimuth 

^Pz 

9  azimuth 

^PsEZ 

5ps 

9  elevation 

^Pe 

9  elevation 

dpz 

9  elevation 

.  9Ps 

dpE 

dpz 

(D.5) 


Equation  ( D.2 ),  the  identity  matrix,  has  no  effect  on  the  H  matrix.  When  equations  ( D.4 )  are  inserted 
into  equation  ( D.5 ),  and  then  multiplied  by  equation  ( D.3  ),  the  final  form  of  the  H  matrix  results. 


9  range 

9  range 

9  range 

9G 

^Ps 

9  azimuth 

^Pe 

9  azimuth 

3pz 

9  azimuth 

9R" 

3ps 

9  elevation 

^Pe 

9  elevation 

dpz 

9  elevation 

.  3Ps 

^Pe 

^Pz 
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Appendix  E  Integrator  Track  Compression  Equations  of  Variation 


This  appendix  develops  the  equations  describing  the  36  elements  of  the  A  matrix  of  the  equation 
below.  These  equations  are  called  the  equations  of  variation.  [7] 


—  0(t,to)  =  A(t)0(t,t„) 


(4.2) 


where 


^(t)“  Vxg|xo(t) 

and  g  comes  from  the  equations  of  motion  of  the  state  vector. 


dt 


X  =  g(X,t) 


The  36  element  form  of  the  A  matrix  is  used  for  the  integrator  track  compression  (fynamics 
model.  Therefore,  the  djynamics  do  not  include  air  drag  and  the  state  vector  does  not  include  the  ballistic 
coefficient.  The  A  matrix  is  defined  as  the  partial  derivative  of  the  equations  of  motion  with  respect  to  the 
state  vector.  Equation  ( 3.8 ),  below,  shows  the  equations  of  motion  for  the  six-element  state  vector. 


X  = 


V, 

31, 

1+o-T 

2 

3  J, 


r 


.2  \ 


1-5- 


H- 


2r' 

3 

1  +  r-f 

2  r 


1-5- 


3-5- 


.2  N 


7J 


(3.8) 


where  x,  y  and  z  are  the  three  components  of  the  position  vector,  ,  Vy  and  are  the  three  components 
of  the  velocity  vector  and  r  is  the  magnitude  of  the  position  vector. 
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Differentiating  equation  ( 3.8 )  with  respect  to  the  state  vector  yields  the  following  36 
components.  Three  of  the  three-%-three  submatrices  are  trivial.  The  remaining  nine  elements  will  now 


be  determined  individually. 


0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

ax 

0 

0 

0 

0 

0 

1 

ax" 

ax. 

ax. 

ax. 

0 

0 

0 

ax, 

ax, 

ax. 

ax, 

ax, 

ax, 

0 

A 

ax, 

ax2 

ax, 

0 

0 

ax, 

ax. 

ax, 

0 

0 

0 

ax, 

ax2 

axj 

(E.l) 


As  in  previous  developments,  aj^lication  of  the  chain  rule  allows  a  simpler  development  as  well  as  easier 
computer  coding. 

dX, 


ax, 


_  X  A  ^  ~  Z 

“  J.5  “  ^3  “  O  ^5 


1  3  ^2 
r"  2  r 


1  +  35 


2  2 
X  z 


-5 


2  ,  2 
X  +Z 


dX, 

ax. 


r>  ^  2  r’  L  r 


(E.2) 


ax, 

ax. 


xz 


xzT  ^  ,  15  105 


3,  T  1^  ,  A-/  Xj 

J+J2  7  l^+z^Z  2 


ax 

ax 


r  -  r'  L  2  ‘  2  r 


1  r’  2  r’  L  r  J 


y'  1  3J2r  y'z'  y'+z' 

- L_32___i._2.^  1  +  35-^ - -5- - 

aX2  r’  r'  2r’L  r"  r" 


ax, 


,,^15  105 

15  +  yZ-^ 


E-2 


ax, 

ax, 


=  3- 


xz  15  J,  xz 


’-’t] 


ax 

ax 


2  r’  L  r  J 

—  =  3—  -—  -  — — r 3-30  — +  35  —1 


Inserting  equations  (  E.2 )  into  equation  ( E.3  )  yields  the  final  form  of  the  A  matrix. 
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Appendix  F  Global  Estimate  Equations  of  Variation 


This  appendix  develops  the  equations  describing  the  49  elements  of  the  A  matrix  of  the  equation 
below.  These  equations  are  called  the  equations  of  variation. 


where 


dt 


a>(t,to)  =  A(t)0(t,to) 


A(t)  =  Vxg 


Xo(t) 


(4.2) 


andg  comes  fiom  the  equations  of  motion  of  the  state  vector, 

|;^x  =  g(x,t) 


The  49  element  form  of  the  A  matrix  is  used  for  the  Global  Estimate  and  Bayes  Filter  dynamics 
models.  Therefore,  the  (fynamics  include  air  drag  and  the  state  vector  includes  the  ballistic  coefficient 
parameter. 


X(7)s 


(5.6) 


The  A  matrix  is  defined  as  the  partial  derivative  of  the  equations  of  motion  with  respect  to  the 
state  vector.  Equation  (  F.  1 ),  below,  shows  the  equations  of  motion  for  the  seven-element  state  vector. 


X  = 


V. 


V, 


Agx 

■^gy  ^iy 


(F.l) 


F-1 


where  V  ,  and  V,  are  the  three  components  of  the  velocity  vector,  Ag  is  the  acceleration  due  to  gravity 

X  y  z 

given  by  equation  ( F.2 )  and  Aj  is  the  acceleration  due  to  air  drag  obtained  from  equation  ( 5.5 )  and 
given  in  component  form  ly  equation  ( F.3  ). 


(F.2) 


where  x,  y  and  z  are  the  three  components  of  the  position  vector  and  r  is  the  magnitude  of  the  position 
vertor. 


0 

0 

0 

-iB*p{v.-(0.y) 

^(V..-co,y)'+(v,  +  ci).x)’ 

■  +  v/ 

-iB»p(v,+ffl,x) 

l/{v,-to,y)’+(v,  +  (ii,x)' 

-|b*pv. 

/( V,  -  tOe  y)^  +  (Vy  +  (0®  x)' 

0 

whereCO®  =0.05883359980154919  [22]  and  the  local  atmospheric  density  is  determined 

using  the  atmospheric  model  described  in  Appendix  A. 
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The  A  matrix  is  formed  from  the  partial  derivative  of  equation  ( F.  1 )  with  respect  to  the  seven- 
component  state  vector. 


ax 

ax 


0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

ax^ 

ax^ 

ax^ 

ax4 

ax^ 

ax^ 

ax^ 

ax, 

axj 

ax3 

ax. 

ax. 

ax, 

ax. 

ax, 

ax, 

ax. 

ax, 

ax. 

ax. 

ax. 

ax, 

ax, 

ax. 

ax. 

ax, 

ax, 

ax, 

ax, 

ax. 

ax. 

ax. 

ax. 

ax. 

aJc, 

ax, 

ax, 

ax, 

ax^ 

ax, 

ax. 

ax^ 

0 

0 

0 

0 

0 

0 

0 

(F.4) 


Of  the  49  elements,  28  are  trivial.  The  remaining  21  will  now  be  determined.  As  in  previous  instances, 
use  of  the  chain  rule  will  simplify  the  notation.  Also,  the  partial  derivatives  of  the  acceleration  terms  Aie 
to  gravity  were  determined  in  Appendix  E  as  equations  ( E.2 ).  Only  the  partial  derivatives  of  the 
acceleration  terms  due  to  air  drag  will  be  developed  here. 

ax 

AA  1  /  \  ap  f/  \2  /  \2  _  _  „ 


axf = "  2  ®  y)  ~  y)' 

1  */  ^  “®p(Vy+CO*x) 

-  y)- - 

[( V,  -  COe  y)'  +  (Vy  +  ©e  x)  +  V,'  ] 

=  *  |^(V,  - ©e  y) [(V,  - ©e  y)'  +  (Vy  -b ©e  x)'  + 


(F.5) 


ax 


ax. 


^1b*p- 


®®[l+(v,-cOe  yf 


ax 


4d 


ax. 


[(v,  -  CO.  yf  +  (v,  +  <0,  x)'  +  V,*  f’ 

=  -jB*  |£(v,  -CO.  y)[{v,  -CO.  yf  +(v,  +co,  x)’  + 
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ax 


4d 


ax. 


=  -^B*p[(v,-0).y)'+(v,+Ci),x)’+V.=]' 


B*p 


r  "I  ^ 

( V,  -  co^  y)'  +  (Vy  +  (0*  x)'  +  V/  ' 


^  =  -Ib*p 
ax,  2“  P 


(v*-co^  y)(Vy +a)^x) 


(V;,  -  cOe  y)'  +  (Vy  +  CO®  x)  + 


ax 


ax,  2 


V,( 

s 

1 

'' 

[(v.-co,y) 

% 

(Vy  +0)eX 

)’+v/ 

Yi 

1 

(V,-C0e  y) 

ax,  ■  2P' 

[(v.  -  (0,  yf  +  (v,  +  00,  x)’  +  V.'  ]  I 

®®K-0>®  y)  (Vy  +  co®x) 


{ V*  -  COe  y)'  +  (Vy  +  (0*  x)'  +  V," 


ax 


5d 


ax,  =  -iB*(v, +“.  ’<)  If  [(v,  -<».  yf  +(v,  +00.  x)’  +  v.’] 


Yi 


+iB*(v,+(».x)- 


®«p(v.-<iJ,y) 


[( V,  -  u>,  y)’ +  (v,  +  (0,  x)*  +  V.’ 

=-^B*  If (v,  +m. x)[(v.  -10, yf  +  (v,  +0), xf  +  V,=]^ 


axj  2  dz 
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ax 


5d  1 


dX,  2 

ax 
lx 


=  -^B*p 


(Vx-coe  y)(Vy+a)ex) 


[(v.  - 0)e  y)'  +  (Vy  +  0)^  x)'  +y"f 


^  =  -|b  *  p  [( V,  - cOe  y)'  +  (Vy  +  0)®  x)'  +  V,'  Y 
_ (Vy +tOe  x) _ 

[(v,-«),yr  +  (v,+(».x)‘+v/]^ 


-iB*p 


8X 


5d  1 


3X.  2 


=  -iB*p 


V.(v,+0).x) 


dX 


5d  1 


1[( V.  -  0),  y)’ +  (v,  +  0),  x)' +  V,* 

^  •  |2  V,  [(V.  -  CO,  y)’  +  (v,  +  00,  x)'  +  V.']^ 

-iB*p(v,+(0,x) 


[(v,-co,yf+(v,  +  oo,x)  +V,’]  J 

(Vy  +  CO®  X) _ 


co.V, 


K 


[(v.-(0.yf+(v,+(0.x)’  +  V.^]' 

=  -^B  ♦  V,  [(v,  -  (0.  y)’  +  (v,  +  (0.  x)’  +  V,= ]' 
+|B‘p(v.-co,y) 


[(v,-co®y)%(Vy+co®x)'+V,^]^ 
^  =  -^B  *  ||  V,  [( V,  -  CO®  y)'  +  (Vy  +  CO®  x)'  +  V/ 


ax 
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3X., 

e 

3 

1 

dX. 

/  \2  /  \2 

4 

- 

(Vx-co®y)  +(Vy+cOex)  +V,  ^ 

V,  (Vy+co®x) 

dX, 

/  \2  t 

D 

. 

K-co^y)  +(v,+o)^x)  +V,^J 

3X„ 

_  JLu  ♦  rt  - 

ax, 

““2®  P 

iVi 

(v,  -  to^  y)'  +  (Vy  +  0)«  x)^  +  V/ 

1/4 

-^B*p  (v,-(0ey)"+(Vy+0)ex)  +V,"J 


dp  3p  dr  dp  X 
dx  ~  dr  dx  dr  r 


dp  dp  dr  dp  y 
dy  ”  dr  dy  ~  dr  r 

dp  _  dp  dr  dp  z 
dz  ~  dr  dz  ~  dr  r 


(F.6) 


The  atmospheric  model  provides  the  gradient  of  the  atmospheric  density  for  use  in  equations  ( F.6 ). 
Inserting  equations  ( E.2 ),  and  ( F.5 )  into  equation  ( F.4  )  yields  the  final  form  of  the  A  matrix. 
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Appendix  G  Orbital  Elements 


The  following  table  provides  the  orbital  elements  for  the  satellites  used  in  simulation  one  and 
simulations  six  and  after.  Note  it  is  given  as  B*,  not  B*/2. 


Table  G-1  Orbital  Elements 


axis 

(DUs) 

ecc 

apogee  perigee 
(DUs)  (DUs) 

inc 

(deg) 

iSB 

IBI 

m^/kg 

1 

1.1431 

0.0649 

1.2173 

1.0689 

66.6670 

337.0477 

149.2720 

111.0444 

0.0871 

2 

1.1176 

0.0207 

1.1408 

1.0944 

71.4925 

33.6281 

24.0132 

299.7730 

0.0978 

3 

1.1189 

0.0260 

1.1480 

1.0898 

71.1933 

87.8515 

239.2868 

14.8873 

0.1107 

4 

1.1284 

0.0422 

1.1760 

1.0807 

17.9057 

84.9663 

288.6909 

157.0855 

0.0326 

5 

1.1023 

0.0207 

1.1251 

1.0795 

52.7924 

239.9506 

13.1372 

186.1254 

0.0554 

6 

1.0882 

0.0009 

1.0892 

1.0872 

5.5404 

37.6862 

69.0430 

254.3590 

0.0337 

7 

1.1030 

0.0374 

1.1443 

1.0617 

72.2489 

41.9233 

119.3513 

128.8672 

0.0365 

8 

1.0933 

0.0134 

1.1079 

1.0786 

4.7114 

265.1745 

77.5812 

201.5614 

0.0191 

9 

1.0744 

0.0398 

1.1171 

1.0316 

58.5559 

81.1987 

224.7762 

133.9783 

0.0605 

10 

1.1655 

0.0231 

1.1925 

1.1386 

39.7510 

148.1973 

72.3493 

14.8886 

0.0613 

11 

1.1305 

0.0028 

1.1337 

1.1273 

10.7826 

66.0109 

51.2450 

88.9736 

0.1042 

12 

1.1594 

0.0438 

1.2101 

1.1086 

62.3430 

218.0545 

30.7839 

247.9762 

0.1302 

13 

1.1165 

0.0516 

1.1740 

1.0589 

2.9282 

282.5804 

39.1168 

8.1137 

0.0096 

14 

1.0627 

0.0243 

1.0885 

1.0369 

33.2378 

240.0768 

346.7636 

65.6831 

0.0739 

15 

1.0948 

0.0327 

1.1306 

1.0590 

31.5112 

206.3120 

204.8270 

166.3494 

0.0378 

16 

1.1828 

0.0771 

1.2740 

1.0916 

13.0270 

353.7673 

100.4942 

314.4835 

0.0052 

17 

1.1267 

0.0061 

1.1336 

1.1198 

8.9446 

17.8261 

84.0984 

153.2555 

0.0474 

18 

1.1697 

0.0145 

1.1866 

1.1528 

43.1829 

171.8751 

20.9602 

169.0522 

0.0548 

19 

1.0947 

0.0230 

1.1199 

1.0695 

13.8640 

14.1571 

95.0683 

123.6986 

0.1150 

20 

1.1209 

0.0667 

1.1957 

1.0461 

50.8408 

266.5706 

147.4709 

289.1914 

0.0815 

21 

1.0681 

0.0063 

1.0749 

1.0614 

27.9660 

183.5601 

311.3900 

27.5289 

0.0142 

22 

1.1650 

0.0143 

1.1816 

1.1483 

58.1834 

287.9121 

4.4881 

285.3072 

0.0339 

23 

1.0869 

0.0315 

1.1211 

1.0526 

27.7908 

275.2574 

71.9475 

57.3137 

0.0762 

24 

1.0769 

0.0304 

1.1096 

1.0442 

64.8925 

262.4626 

214.3987 

319.7559 

0.0203 

25 

1.1694 

0.0884 

1.2728 

1.0661 

6.1776 

102.7449 

238.3270 

121.0557 

0.1074 

26 

1.0541 

0.0112 

1.0659 

1.0423 

16.7671 

69.6316 

227.7882 

84.9530 

0.0274 

27 

1.0708 

0.0291 

1.1020 

1.0397 

52.3488 

28.1130 

204.8310 

359.6201 

0.0012 

28 

1.1139 

0.0084 

1.1233 

1.1046 

49.0407 

222.6593 

189.3927 

333.3271 

0.0698 

29 

1.0553 

0.0165 

1.0727 

1.0379 

4.8674 

327.9385 

0.3699 

216.3699 

0.1252 

30 

1.0844 

0.0357 

1.1232 

1.0457 

33.3252 

320.9543 

184.2456 

277.3348 

0.0372 

31 

1.2088 

0.0181 

1.2306 

1.1869 

50.9498 

301.9593 

197.2199 

163.3539 

0.0789 

32 

1.1266 

0.0175 

1.1462 

1.1069 

11.8503 

40.0213 

108.4510 

277.6716 

0.0002 

33 

1.1065 

0.0021 

1.1089 

1.1042 

1.9064 

281.0339 

268.8070 

102.5428 

0.0050 

34 

1.1203 

0.0096 

1.1310 

1.1095 

12.8023 

297.6429 

281.1327 

356.7522 

0.1513 

35 

1.1373 

0.0056 

1.1436 

1.1310 

7.8831 

2.4541 

20.7305 

242.0845 

0.0312 

36 

1.1360 

0.0551 

1.1986 

1.0735 

18.4088 

189.6166 

344.1614 

95.7255 

0.1098 
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37 

1.1555 

0.0801 

1.2479 

1.0630 

26.7252 

1.8692 

228.3934 

138.8889 

0.0736 

38 

1.0863 

0.0478 

1.1382 

1.0344 

7.3939 

148.0673 

274.4587 

36.7271 

0.0954 

39 

1.0781 

0.0211 

1.1009 

1.0554 

30.0380 

314.0126 

264.2481 

1.7241 

0.0323 

40 

1.0992 

0.0217 

1.1231 

1.0753 

65.3707 

31.9341 

330.4760 

86.7126 

0.1797 

41 

1.1745 

0.0087 

1.1847 

1.1642 

8.7431 

269.2701 

138.8004 

72.5979 

0.0946 

42 

1.1206 

0.0658 

1.1944 

1.0469 

17.7696 

164.5883 

2.1040 

43.9168 

0.0466 

43 

1.1186 

0.0465 

1.1706 

1.0666 

28.7028 

75.3733 

183.8061 

304.7195 

0.0772 

44 

1.0808 

0.0220 

1.1046 

1.0570 

37.5336 

26.2564 

223.7465 

107.4421 

0.0924 

45 

1.0904 

0.0184 

1.1105 

1.0704 

53.9295 

4.1245 

12.5513 

342.3633 

0.0102 

46 

1.1576 

0.0161 

1.1762 

1.1389 

59.6470 

131.5987 

202.4999 

153.0539 

0.0972 

47 

1.0882 

0.0254 

1.1158 

1.0605 

24.9452 

278.1446 

117.9650 

168.0659 

0.1082 

48 

1.0660 

0.0132 

1.0800 

1.0520 

45.4555 

245.8460 

190.7447 

76.6971 

0.0755 

49 

1.1678 

0.0164 

1.1870 

1.1486 

48.6396 

27.0271 

345.1981 

99.4496 

0.0828 

50 

1.0816 

0.0232 

1.1067 

1.0566 

8.0752 

9.2708 

29.5075 

195.4267 

0.0490 

51 

1.0739 

0.0088 

1.0833 

1.0645 

8.2498 

319.9620 

165.1917 

329.4712 

0.1095 

52 

1.1033 

0.0287 

1.1350 

1.0716 

20.9109 

320.3864 

132.8557 

39.8585 

0.1393 

53 

1.0633 

0.0101 

1.0740 

1.0526 

55.6062 

164.4913 

1.1809 

121.1371 

0.0917 

54 

1.0537 

0.0193 

1.0740 

1.0334 

17.2056 

118.1681 

250.3306 

307.3486 

0.1282 

55 

1.0876 

0.0425 

1.1338 

1.0413 

46.6057 

112.6308 

81.9692 

337.5289 

0.0069 

56 

1.1219 

0.0335 

1.1594 

1.0843 

4.8253 

345.4804 

46.4249 

337.4609 

0.0347 

57 

1.0431 
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